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Exactly soluble models in the theory of electromagnetic propagation and scattering are essentially 
restricted to horizontally stratified or spherically symmetric geometries, with results also available 
for certain waveguide geometries. However, there are a number of new problems in remote sensing 
and classification of buried compact metallic targets that require a wider class of solutions that, if 
not exact, at least support rapid numerical evaluation. Here, the exact Chandrasekhar theory of the 
electrostatics of heterogeneously charged ellipsoids is used to develop a "mean field" perturbation 
theory of low frequency electrodynamics of highly conducting ellipsoidal targets, in insulating or 
weakly conducting backgrounds. The theory is based formally on an expansion in the parameter rjs = 
Ls/Ss{ijj), where La is the characteristic linear size of the scatterer and Ss{ljj) is the electromagnetic 
skin depth. The theory is then extended to a numerically efficient description of the intermediate- 
to-late-time dynamics following an excitation pulse. As verified via comparisons with experimental 
data taken using artificial spheroidal targets, when combined with a previously developed theory 
of the high frequency, early-time regime, these results serve to cover the entire dynamic range 
encountered in typical measurements. 



I. INTRODUCTION 

There are a number of longstanding economic and hu- 
manitarian problems, such as clearance of unexploded 
ordnance (UXO) from old practice ranges, that require 
remote identification of buried metallic objects. The 
most difBcult technological issue is not the detection of 
such targets, but rather the ability to distinguish between 
them and harmless clutter items, such as various sized 
pieces of exploded ordnance. Since clutter tends to ex- 
ist at much higher density, even modest discrimination 
ability leads to huge reductions in the economic cost of 
remediating such sites. 



A. Electromagnetic inverse problems 

Formally, a successful solution to the electromagnetic 
(EM) discrimination problem is a theory or algorithm 
that allows derivation of accurate bounds on physical 
properties of the target scatterer (its position, shape, 
orientation, physical composition, etc.) from measure- 
ments of the scattered field using a well characterized 
experimental apparatus (with known transmitter and re- 
ceiver coils, transmitted waveform, and so on). Solu- 
tion of this inverse problem first requires the ability to 
generate high-fidelity candidate solutions to the forward 
problem, namely accurate computations of the scattered 
field from a known target in a known subsurface envi- 
ronment. The general solution to the forward problem 
requires full three dimensional numerical solutions to the 
Maxwell equations, a difficult and time consuming com- 
putational problem. To reduce the computational bur- 
den, it is extremely important to obtain analytic solu- 
tions to as broad an array of exactly soluble model prob- 
lems as possible. These solutions may then either be used 



as crude models of the target, or as the basis of a pertur- 
bation scheme for accurate modeling of "nearby" target 
geometries. 

The only compact targets for which a full analytic so- 
lution at any frequency may be derived are those with 
spherical symmetry [1]. These are rather poor approx- 
imations to UXO, which tend to more resemble finite, 
rounded cylinders with roughly 4:1 aspect ratio. The ap- 
proach pursued here is to take advantage of the fact that 
electrostatic solutions exist for a much broader array of 
target geometries, and that these solutions can then be 
used as the basis for a controlled perturbation theory, 
valid at low frequencies. The small parameter in the the- 
ory, = Ls/Ss, is the ratio of the electromagnetic skin 
depth Sg{uj) to the linear target size L^. The theory is 
dubbed the "mean field approach," since the smallness of 
77s means that the expansion is highly nonlocal in space, 
with the currents and fields at any given point in the 
target being sensitive to their values throughout the tar- 
get. Although formally valid only for small r^s, we will 
see that the theory may be extended to higher frequen- 
cies, even where rjg is significantly larger than unity, if 
one generates a sufficient number of terms in the series 
[2]. 

The basic zeroth order theory requires one to solve 
for the electrostatic field generated by the target in a 
sequence of background fields of increasing complexity. 
The perturbation theory is developed formally for a gen- 
eral target shape, but even this sequence of simpler elec- 
trostatic problems generally requires a numerical solu- 
tion. However, for the case of ellipsoids, such solutions 
may be computed analytically via an elegant approach 
developed by Chandrasekhar [3, 4]. Since many targets 
of interest may be modeled quite accurately by ellipsoidal 
or spheroidal shapes, the results of this paper have an im- 
mediately relevant application. The theory will mainly 
be illustrated for the case in which both background and 
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scatterer are nonmagnetic (i.e., the permeability is a 
uniform constant), but the extension to permeable tar- 
gets will be described as well. 

When treating buried targets, the electrodynamics of 
the soil is also potentially important. We will assume 
that the ground is insulating or sufficiently weakly con- 
ducting that its response may be treated as quasistatic 
in the frequency range of interest (say, 100 kHz or less). 
Specifically, the background EM penetration depth (typ- 
ically tens of meters or more at these frequencies) should 
be large compared to the measurement domain (typically 
on the scale of 1 m). 

Even with the quasistatic assumption, the total elec- 
tric field has a significant, unpredictable variability due to 
strong variation in the dielectric function due to varying 
soil type and inclusions, surface vegetation, air-ground 
interface, etc. It transpires, however, that an induction 
loop (EMI) measurement (as opposed, say, to a linear 
antenna measurement) is effectively sensitive only to the 
"magnetic part" (curl component) of the electric field, 
and that the latter is insensitive to the ground (par- 
tially explaining the ubiquity of such measurements in 
geophysics), so long as it is nonmagnetic [5]. In a happy 
confluence of theory and experiment, the perturbation 
theory is most efficiently formulated to isolate and com- 
pute precisely this part of the field. 



B. Time-domain measurements 

A common experimental probe for metallic targets is 
the time-domain electromagnetic (TDEM) measurement, 
in which one detects the inductive response of a target 
following termination of a transmitted pulse. The pulse 
generates a characteristic pattern of currents in the tar- 
get, and the subsequent dynamics of, say, the electric field 
may be written as a superposition of EM eigenmodes: 

oo 

E(x,t) -^A.e(")(x)e-^"*, (1.1) 

n=l 

in which e^") is the mode shape, A„ the decay rate, and 
An the excitation amplitude. This is analogous to the re- 
sponse of a drumhead following a strike. However, rather 
than corresponding to a set of characteristic eigenfre- 
quencies (with, perhaps, some weak damping), the dissi- 
pative/ohmic dynamics leads here to modes that exhibit 
a pure exponential decay in time. The data will typically 
consist of the voltage measured in a receiver coil, 

oo 

V{t) = Y,Vne-'-\ (1.2) 
ji=i 

which is a superposition of the same exponential decays. 
The decay rates and mode shapes are, respectively, eigen- 
values and eigenfunctions of the Maxwell equations at 
imaginary frequency a;„ = — iA„, and may therefore be 



accessed through the perturbation technique developed 
here. 

It transpires that in the high-contrast limit there are, 
in fact, two classes of excitation with widely separated 
decay rates. One set corresponds to electric polarization 
of the target, the primary example being an excitation 
in which a uniform electric field is suddenly switched off. 
These "electric modes" relax essentially instantaneously 
by direct equilibration of the charges induced on the tar- 
get surface, and hence have very large, perturbatively 
inaccessible A„. On the other hand, the rapid relaxation 
implies that their contribution to the signal (1.2) disap- 
pears almost immediately. 

The second, more interesting set of modes, corresponds 
to magnetic polarization of the target, the primary ex- 
ample being an excitation in which a uniform magnetic 
field is suddenly switched off. In this case circulating cur- 
rents are generated. The flows are essentially tangential 
at the target boundary, produce no charge polarization, 
and therefore take much longer to relax. Being intrinsic 
to the geometry of the target, the current patterns associ- 
ated with these "magnetic modes" must vary on the scale 
Ls (or even on much smaller scales, as the decay rate in- 
creases and the mode shape becomes more spatially com- 
plex), and hence must lie in the regime 77s(— *A„) = 0(1). 
Once again, however, these may be accessed via a suffi- 
ciently high-order expansion [6]. Moreover, the exponen- 
tial decay implies that as time progresses fewer and fewer 
of even these modes contribute to the signal, so that such 
a theory, which accurately computes only a finite set of 
the slowest decaying modes, would provide quantitative 
predictions on intermediate-to-late-time scales. On the 
other hand, in the opposite, early-time regime when a 
very large number of modes simultaneously contribute 
(but not so early that any electric modes still survive), 
a complementary "surface mode" theory, based on the 
diffusion of the initial screening current inward from the 
target surface, may be developed [7, 8]. It will be seen 
that these two regimes significantly overlap for a suffi- 
ciently high order mean field expansion, enabling a quan- 
titative prediction of the signal over the full time-domain 
dynamic range. 



C. Outline 

The outline of the remainder of this paper is as follows. 
In Sec. II the basic content of the method is introduced, 
showing that it reduces to the evaluation of certain inte- 
grals of the Coulomb interaction over the volume of the 
scatterer. In Sec. Ill the Chandrasekhar theory of ellip- 
soidal electrostatics [3] is reviewed, showing that it too 
reduces to this same class of integrals. In Sec. 1\ this con- 
nection is used to derive explicit expressions for the scat- 
tering integrals. In Sees. A' and Yl these are evaluated 
explicitly for the case of solid homogeneous ellipsoids. 
In Sec. \ II the solution to the frequency-domain scat- 
tering problem in a known background field is obtained 
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first formally, and then explicitly to 0{ri1). Electric and 
magnetic excitations are identified, the latter being the 
low frequency precursors to the freely decaying magnetic 
modes. In Sec. Mil, the time-domain response is dis- 
cussed. Successful comparisons of the theory, extended 
numerically to significantly higher order in ?7s, to real 
data from artificial spheroidal targets are demonstrated. 

Finally, in Sec. we conclude by describing general- 
izations of the theory (whose detailed applications will be 
left to future work) . In Sec. IX A we consider permeable 
targets, /i ^ ^b. In Sec. IX B we consider simplifications 
in the high contrast limit /i//J.b ^ 1, relevant to ferrous 
targets where n/fib — O(IO^). In Sec. IXC we consider 
the computation of the freely decaying modes for per- 
meable targets. In Sec. IX D we consider more realistic 
target geometries, including hollow targets and multiple 
targets. Finally, in Sec. IX E we consider the effects of 
background permeability variations. 



II. MEAN FIELD APPROACH TO LOW 
FREQUENCY, HIGH CONTRAST SCATTERING 

The mean field approach is based on the Green func- 
tion formulation of electromagnetic scattering. Let 
eti(x, cli) and e(x, oj) be the space- and frequency- 
dependent dielectric constant for background and back- 
ground plus target, respectively. At low frequencies these 
take the form (in Gaussian units, which are used through- 
out unless explicitly stated otherwise), 

£6 = 4 + e = e'+ , 2.1 

where e[,(x), e'(x) are static dielectric constants and 
crh(x), cr(x) are DC conductivities. All four quantities 
are real and frequency independent. Let k = uj/c he the 
vacuum wavenumber, and define 



Kf, — ebiJ-bk , K = e/ifc , Q = K 



(2.2) 



For most of this paper, it will be assumed that the per- 
meability /i = /ift is a uniform constant, the same for 
both background and target. Generalization to inhomo- 
geneous /i will be discussed in Sec. IX. In addition, at 
the low frequencies of interest, e'{,,e are negligible com- 
pared to the conductivity contributions and will usually 
be dropped [9]. The high contrast assumption corre- 
sponds to = \e/€b\ ^ (j/(7b ^ 1- Typical metalhc 
conductivities are in the range a ^ 10^ S/m, while typi- 
cal ground conductivities are in the range ab ^ 0.1 S/m, 
so the ratio a/ab ^ 10^ is indeed extremely large. 



A. Green function formulation 

The background and full electric fields (with uniform 
fi) satisfy the frequency domain wave equations 



V X V X Efc - 
V X V X E 



- K^E 



S 



in which S(x) = (47ri/ifc/c)js(x) is proportional to the 
source current distribution jg. The latter is generated by 
a transmitter coil that is assumed to be external to the 
target region. The associated source charge distribution 
is PS — {iuj)~^W js- The magnetic fields follow from the 
relation ifcB — ikpR = V x E. 

Let G(x, x') be the 3x3 tensor Green function for the 
background medium, satisfying, 

V X V X G(x, x') - /tf,(x)2G(x, x') (5(x - x')!. (2.4) 

G is symmetric, but self adjoint only if is real. For a 
uniform, homogeneous background medium one obtains 
the exact solution 

G(x, x') = (^1 + -i VV^ .9(x, x'), (2.5) 



with scalar Green function 



5(x,x') 



giK[,|x-x'| 

47r|x - x'l 



(2.6) 



satisfying the Helmholtz equation, — (V^ + K^)g = (5(x — 
x'). The background electric field may then be expressed 
in the form 



Efc(x) = / <i3x'G(x,x') • S(x'), 



(2.7) 



The operator acting on g in ( >) essentially enforces 
transverse polarization via the divergence condition V • 
(ejjE;,) = ^T^Ps- Note that quite generally, as indicated by 
(2.5) and (2.6), G has a |x— x'j^'' divergence at small sep- 
aration, hence there is a potential logarithmic singularity 
in ( ) . The integral over the anisotropic dipole-like an- 
gular dependence regularizes this singularity, but careful 
limiting procedures must still be used when dealing with 
Green function integrals of this type. 

By taking the difference of the two lines of ( ), the 
source S drops out and one obtains 



V X V X (E - Eb) - K^(E - Eb) = QE. 



(2., 



Applying the Green function to both sides, the full field 
satisfies the integral equation, 



E(x) = Eb(x) 
= Eb(x). 



d3a;'Q(x')G(x,x') •E(x') 

1 + ^VV ) • / d^a;'.g(x,x') 

X Q(x')E(x'), (2.9) 



(2.3) 



where second line is valid for the case of a homogeneous 
background. The ability to factor the double gradient 
outside of the integral renders the latter explicitly con- 
vergent in this case. The utility of this formulation is 
that Q(x') vanishes outside of the scatterer volume, Vs- 
Hence if x is restricted to 14, (2.!)) becomes a closed equa- 
tion for the field internal to the scatterer. When x lies 
outside of I^, the external field then follows by simple 
integration of the internal field. 



4 



B. High contrast limit 

As alluded to earlier, even though £{, is assumed very 
small, it can still vary by many orders of magnitude (e.g., 
between weakly conducting ground and insulating air), 
generating highly variable contributions to G and to the 
background electric field. However, we will now see that 
this component does not contribute to a magnetic field 
or EMI measurement (though it will contribute, e.g., to 
a linear antenna electric field measurement), and can be 
removed from the computation at the outset. This is a 
key result because the non-inductive part of the field, due 
to this intrinsic variability, is essentially unpredictable. 

The fact that the — limit of ( ) is singular 
is evident from the fact that QE is not divergence free 
(at minimum, there is a delta-function contribution due 
to the target boundary discontinuity); this is also evi- 
dent from the term on the right hand side of (2.9). 
Therefore (2.8) has no solution if the kI term on the left 
is simply dropped. To account for this, decompose E and 
Eb into divergence free and curl free parts, 

Eb = ikAb - V$f, 
E = ifcA-V$. (2.10) 

Here the vector potentials generate the magnetic field via 
V X A = B, V X Ab = Bb and, for later convenience, we 
choose the Coulomb gauge 

V-A = 0, V-Ab = 0. (2.11) 

It will be seen that a consistent Kb-independent solution 
for A, Ab exists when Kb — >■ 0, while continue to 

depend strongly on Kb in this limit. However, the inte- 
gral of a gradient around any closed curve vanishes, and 
this variability then indeed disappears in an EMI mea- 
surement. 



1. Background field 

Consider first the background field. Using (2.11), the 
first line of (2.3) may be written the form [9] 

_V2Ab = ^(js+abEb), (2.12) 
c 

with formal solution 

Ab(x) = ^ / ,3^ Js(xO + ab(x-)Eb(xO ^ 

This result is equivalent to the Biot-Savart law for 
the magnetic field arising from the combination of 
source/transmitter current js and the induced back- 
ground currents ab^b- For an insulating or weakly con- 
ducting background the latter term is expected to be neg- 
ligible, and the formal limit Kb ^ yields 

Ab(x) = ^ / dVi^, (2.14) 
cJ |x-x'| 



which is indeed entirely independent of the background 
details. 

To establish conditions for consistency of this conclu- 
sion, an equation for $b is obtained by taking the diver- 
gence of both sides of ( ) : 

V • (CTbV$b) = -iujps + ik{V(7b) ■ Ab. (2.15) 

The first term on the right generates the usual static 
Coulomb field (distorted by the nonuniform eb). If ps 
is nonzero, then (TbEb could indeed be of the same or- 
der as is, and (2.14) is invalid (the CbEb ~ — crbV$b 
term contains a potentially large quasistatic Coulomb 
contribution). However if, as is typical, the transmit- 
ter is purely inductive, i.e., does not generate any free 
charges, ps — iwV • js = 0, then this term is absent, 
(2.14) is valid, and ( ) may put in the form 

— V-(abV$b) = -(Vab)-Ab. (2.16) 

CTb Cb 

Both sides of this (generalized Poisson) equation for $b 
depend explicitly on the "shape" of Kb(x), but not its 
overall magnitude. Substituting (2.14) into its right hand 
side, the formal Kb ^ limit then produces a finite value 
of $b, but with, as alluded to earlier, a complicated spa- 
tial dependence that depends on the detailed geometry of 
the ground conductivity. The physical interpretation of 
this result is that the oscillating magnetic field generated 
by (2.14), in combination with the nonuniform conduc- 
tivity, induces a small background charge density, which 
then produces a finite quasistatic electric field. However, 
the contribution — (TbV<I>b to the background current for- 
mally vanishes when Kb — >■ 0, producing (2.14). 

Note that, comparing (2.ri) and (2.6), this limit also 
implies a near-field approximation in which propagating 
wave effects are neglected. Quantitatively, this requires 
that |Kb|i? <C 1 where R is the length scale of the mea- 
surement domain (e.g., transmitter-target separation). 
This is equivalent to R/^ <^ 1, where f = c/y^2nfiabUJ 
is a characteristic background skin depth. In MKS units 
one obtains 

e _ 5 [PoV^^ / O.l S/m y/' /lOkHzy/' 

which confirms that in the parameter domain of interest 
the approximations we have used will certainly be valid 
for R on the scale of a few meters or less. 

We emphasize again that the contribution of Eb to 
an induction measurement is only through Ab, which is 
insensitive to the background conductivity (and, more 
generally, to [9]) as claimed. 

2. Full field 

We now proceed in a similar fashion to obtain a con- 
venient formulation of equation ( ' '') for the full electric 
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field in the high contrast hmit. Analogous to (2,12), sub- 
stituting (2,11)) into (2.8) yields 



A(x) - Ab(x) = 



,3 ,a(x')E(x')-a,(x')E,(x') 



c J |x — x'l 

(2.18) 

Consistently, from (2..'>), the divergence of the right hand 
side vanishes. One expects the scattered and background 
fields to be of the same order in the measurement region, 
and so in the limit ki,/k ^ one obtains, analogous to 
(2.14), 

ifc[A(x) - Ab(x)] = E(x)-Eb(x)+V[$(x)-$b(x)] 



ikfj. 



,a(x')E(x') 



(2.19) 



Since a = at outside Vs, correct to the same order, we 
have also dropped the exterior contribution from the crE 
term and restricted the integral to Vg- 

A generalized Poisson equation for <I> is obtained from 
the divergence of the second line of (2.;>): 



V-(aV$)^-(Va)-A 



(2.20) 



while entirely avoiding the computation of $,$6. The 
basic strategy is to first solve for the full electric field 
Eint inside Vs , and then insert this into the left hand side 
of ( ) to compute A outside Vg. A measurement of 
the magnetic field B = V x A, or of the induced voltage 



V = (h ^■d\ = ik j) A-d\ (2.21) 

'Cr *^Cr 



in a receiver loop Cr, are then independent of <&. 

The key to solving for Ei„t is to project (2. :i!) onto 
the appropriate vector function space, namely the space 
of functions E restricted to the domain Vg and satisfying 



V • (crE) =0, xeVs 
n-E = 0, xedVs. 



(2.22) 



A purely inductive transmitter, ps = 0, has again been 
assumed. The solution to this equation again depends 
on the detailed geometry of cr = (Tb outside Vg. For a ho- 
mogeneous target and background, the right hand side is 
supported at the discontinuity of cr at the target bound- 
ary 



C. Solution strategy 

The remainder of this paper will be focused on solv- 
ing ( ' '") for A [with A;, already determined by (2.14)], 



The first condition is equivalent, in the high contrast 
limit, to V • (k^E) = 0, while the second, in the same 
limit, follows from the usual EM boundary condition, 
namely continuity of en • E across the boundary dVg, 
where n is the local surface normal. The latter implies 
that n • Eint = (ef,/e)n • Ecxt — >■ 0. Physically this fol- 
lows from the fact that since background currents are 
much smaller than those in the target, the latter must 
fiow parallel to the boundary in order to maintain charge 
conservation. 

Next, the identity 



/ d3xCT(x)E(x) • V/(x) = - / d^xf{ii)y ■ [c7(x)E(x)] + / dAcr(x)/(x)n(x) • E(x) 



= 0, 



(2.23) 



which follows by applying Green's theorem, and is valid 
for any scalar function /, shows that the function space 
(2.22) is orthogonal to the space of gradients, in the sense 
of the above inner product with kernel cr. Denoting the 
orthogonal projection onto the space (2.22) by Va, equa- 
tion (2.19), when restricted to the scatterer volume, takes 
the form 

^/ N ^ ^ . X ika-- f ,^ ,cr(x')E(x') 
E(x) - VM^) - —T^a / rfV ^ \.' . (2.24) 
c Jv, |x-x'| 

This is a closed equation for the internal electric field [re- 
stricted to the space (2.22)], and is explicitly independent 
of $b and of the form of $ outside Vg ■ 



1. Basis function expansion 

A more explicit form of this equation, which is the 
foundation for the mean field perturbation scheme to be 
developed, is obtained by expanding 

a(x)E(x) = ^aMZM(x) (2.25) 

M 

in terms of a complete (though not necessarily orthogo- 
nal) set of basis functions Zm consistent with (2.22): 

l^'y^'^ytal (2.26) 
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By inserting this expansion into (2.24), and taking the 
inner product on the left with Z^, one obtains a formal 
matrix equation for the coefficients 



M 

in which 



M 



n - I .3 ZlW*-Zm(x) 



are basis function inner products, 



(2.27) 



(2.28) 



2. Basis functions for ellipsoids 

For general targets, explicit forms for the Zm (or, 
equivalently, for the projection operator Va) may be hard 
to come by, but for ellipsoids they may be constructed 

directly from known forms for the sphere. Specifically, 

( s) 

if Z\/{x) are basis functions obeying (2.26) on the unit 
sphere, then 



Z'y\x.) = Qa ea Z[^^^ (xi/ai,X2/a2,X3/a3), (2.32) 



XT / ^3 / j3 /Zl(x)* • Za/(x') 



(2.29) 



represents the projection of the Coulomb integral onto 
the basis function, and 



ab,L 



[ d3^Zi(x)* •Efc(x) 
Jv„ 



(2.30) 



are basis functions for the ellipsoid with principal radii 
a = (ai, 02, 03), where is the unit vector along princi- 
ple axis a. The unit surface normal is given by 



n(x) = go^aK ^^QVs, 



(2.33) 



represents a similar projection of the background field. 
The V(<i> - <&{,) term drops out via (2.23). The formal 
inverse 



a=(0-^H 



a;, 



(2.31) 



where the boundary is defined by 



(2.34) 



determines the internal field expansion coefficients in 
terms of those for the background field. Since a is real, 
the matrices O, H are both self-adjoint. 



A convenient complete set of basis functions for the 
sphere may be constructed from the vector spherical har- 
monics [1] 



1 



(2.35) 



in which L = — zxx V is the angular momentum operator 
and cirn = •\/(/ — m){l + to + 1). The key properties are 
that the vector harmonics are tangent everywhere to the 
sphere surface, and divergence free. 





0, 



and the functions Xi„j, x x X/™, I > 1, —I < m < I, form 
a complete basis for the tangent vector fields, obeying 
the orthonormality relations 



Su'S„ 



df^Xr„ • (x X Xp,„0 = 0. 



(2.37) 



Furthermore, since x^Yim are polynomials of degree I in 
xi,X2,x^, so are the components of x^~Kim. Similarly, 
x'"''^xx Xi„i is a vector polynomial of degree l+l. Finally, 
note the identity 



(2.36) ^ X [/(2;)X/, 



J-) Jfm H X X A;^ 



(2.38) 

which shows that the choice f{x) — (1 — x^)a;' produces a 
divergence free polynomial of degree I -\- 1 that is tangent 
everywhere to the unit sphere surface, x — 1. 

From the completeness property it follows that a com- 
plete set of basis functions Z^^^^ for the sphere, in the 
form of divergence free vector polynomials that are tan- 
gent at the sphere surface, and in which M = {l,m,p,i) 
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is now a composite index, may be taken as 

Irnp 



Z|™p = Vx [(l-:r2)x'+^PX,„] (2.39) 
= (1 - x^) V X [a;'+2PX,„] - 2x'+'P+'i x X^^ 

for p = 0, 1, 2, . . .. The i = 1 basis functions are of degree 
n = I + 2p > 1, while those for i — 2 are of degree 
n ~ I + 2p + 1 > 2. At a given degree n > 1, there 
are {n + l){n + 2)/2 polynomials of the first type for n 
odd, and n(n + 3)/2 for n even; and there are n{n + 1) /2 
polynomials of the second type for n even, and {n—l){n+ 
2)/2 for n odd. Including both types, there are P{n) — 
n{n + 2) basis functions of degree n for both even and 
odd n. 

The transformation (2.32) produces the correspond- 
ing basis functions zj^^ for the ellipsoid, which are 
then polynomials of identical degree, but with coefficients 
rescaled by appropriate powers of the Oq,. 



D. Time domain eigenvalue equation 

Equation (2.31) describes the target response to an 
external source at fixed frequency. In time-domain mea- 
surements one is instead interested in the free evolution 
(1.1) of the system following pulse termination. Although 
the excitation coefficients A„ will depend on the details 
of the pulse (and their computation will be dealt with in 
Sec. ' )> the mode shapes e^") and decay rates A„ do 
not. The eigenvalue equation these satisfy corresponds 
to the second of equations (2.3) with the replacement 
w = —iX and with vanishing source, S = 0: 



V X V X e(") - K2(-iA„)e("' = 0. 



(2.40) 



Correspondingly, these modes satisfy the homogeneous 
form of the integral equation (2.1 'I) or (2.24) in the ab- 
sence of the background field: 



e(")(x) 



,(7(x')e(")(x') 



(2.41) 



The corresponding basis function expansion of the modes 
a(x)e(")(x)=5]4';)ZM(x) (2.42) 



M 



then produces the self-adjoint generalized eigenvalue 
equation 



(2.43) 



The expansion coefficients may be normalized to obey 
the orthonormality condition 



a('»)tOa(") 



(2.44) 



It is precisely this form of the equations that will be ana- 
lyzed in Sec. VHI. The normalization condition is equiv- 
alent to 



(x) 



(x) - d.„ 



(2.45) 



which follows directly from self-adjointness of the dou- 
ble curl operator in (2.4U). In the high contrast limit, 
the integral may be restricted to Vs, with errors of order 

CTfc/cr. 



E. Role of the Chandrasekhar theory 



The vector harmonics allow one to diagonalize the sys- 
tem (2.>i L) for the sphere, but no such simplification oc- 
curs for more general target geometries. One is therefore 
forced to develop approximate solutions based on compu- 
tation of the array elements (2.28)-(2.3(J) for a truncated 
set of basis functions. 

The key observation, however, is that for homogeneous 
ellipsoids the polynomial character of the basis functions 
allows one to perform the Coulomb integrals ( ■ ) ana- 
lytically. Specifically, one requires integrals of the form 



/ d' 



I mi / ra-i i 



(2.46) 



in which the domain of the x' integral is restricted to the 
ellipsoid interior, but x may be either inside or outside. 
The problem of computing 2?m(x) is isomorphic to that of 
computing the electrostatic potential due to an ellipsoid 
with static charge density p(x) = x'^'^x^^x^^. In his 
book [3], Chandrasekhar presents an elegant formalism 
for computing potentials of precisely this type. Moreover, 
for X £ Vs, I?m(x) is also a polynomial (of degree mi -I- 
rn2 + rn^ + 2). The x-integral in (2.29) therefore has a 
pure polynomial integrand and is trivial to perform. In 
Sec. Ill an overview of his method is presented. Following 
that, the results of these evaluations will be applied to 
the solution of the scattering problem. 

In numerical applications, the matrix equation ( ) 
will be truncated at finite order. As the index M in- 
creases [more specifically, as the indices l,p in (2.39) in- 
crease] , the spatial complexity of the basis functions in- 
crease, so this truncation works best for smoother (typ- 
ically, lower frequency) field distributions. Correspond- 
ingly, mode complexity increases as the decay rate A„ 
increases, and the truncated eigenvalue equation (2.43) 
will be accurate only for a finite set of more slowly de- 
caying modes. 
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III. THE CHANDRASEKHAR THEORY OF 
ELLIPSOIDAL ELECTROSTATICS 

We will now develop a theory for the analytic calcula- 
tion of Coulomb integrals of the form 

for a certain class of charged densities p, which include 
the monomial forms (2.45), for the case in which Vs is an 
ellipsoid, defined by its principal axes a = (01,02,03). 

A. Homoeoids 

For any given ellipsoid one may define a family of sim- 
ilar concentric ellipsoids by 

E 4 = ^' (3-2) 

with /X > 0. Equation (3.2) clearly defines an ellipsoid 
with principal axes ^/Jlai, y/Jla2, y/Jia^- A homoeoid is 
defined to be a shell bounded two similar concentric el- 
lipsoids, i.e., the set of points x such that 

m<E#^^2, (3.3) 

a " 

for some < A*i < /^2- An homogeneous homoeoid is one 
that has a uniform charge density over its interior. An 
infinitesimal homoeoid is the 2D surface that emerges 
when /i2 — /^i — > 0. A homogeneous infinitesimal ho- 
moeoid will be called a hi-homoeoid. Initially a special 
kind of heterogeneous homoeoid will be considered, in 
which surfaces of constant charge density are concentric, 
similar hi-homeoids: p = p{p) only. These will be called 
s-homoeoids (or s- ellipsoids if there is no inner bounding 
surface). In particular, p in (3.1) will be taken to be a 
function of p, alone. Since one may view a s-homoeoid as 
a superposition of hi-homoeoids, the potential due to a 
s-homoeoid may be written in the form 

</)(x)= f' p{p)cl>{p-^)dp (3.4) 
J 111 

in which (/)(/!; x)(i/i is the potential due to a hi-homoeoid 
of uniform unit charge density bounded by p and p + 



dp. Moreover, by simple homogeneity of the Coulomb 
integral one obtains 

J iM-^/^x-y'l 1^ ^al) 

= 0(1; x/V^), (3.5) 
where the change of variable x' = ,/p.y' has been used. 
Thus, all results for s-homocoids may be obtained by 
considering only the potential due to a hi-homoeoid with 
parameter p = \. 



B. Potential in the interior of a hi-homoeoid 

Consider first the potential interior to the hi- 
homoeoid. Note that (by the usual rules for manipulat- 
ing arguments of (5-functions) the surface charge density 
implied by ( ) is nonuniform, even though the origi- 
nal volume charge density is uniform. Using spherical 
coordinates with origin at the observation point x, one 
obtains 

0(l;x)- j dnJ\drS[f{r)]^ J dn^j^^^, (3.6) 
with 

m ^ ""2'"°^' - 1' (3-7) 

where n{n) is the unit vector in the direction fl, r{^l) is 
the positive solution to /(r) = 0, and the 1/|/'| factor 
follows from the rules for integrating delta-functions with 
nontrivial arguments. Recalling that ^a/o-a = ^/R^ i 
where R{^) is the radius of the shell in direction n mea- 
sured from the natural origin at the ellipsoid center, one 
obtains a quadratic equation for r(f2): 

= r^-)-2i?^.;^^-Ki?^fE^-l), (3.8) 

with solutions 



± 



i?4 



(3.9) 



Since X^a^a/'^a < 1 (x is interior to the shell), these solutions exist for any fi. One may now write /(r) 
(r — r+){r — r^)/R^, and hence f'{r+) — —f'{r-) = {r^ — r^)/R^. The positive solution is r+, while r_(f2) 
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— r_|_(— f2) corresponds to the point on the eUipsoid in the direction — fi. Averaging these two contributions and using 
R{fl) = R{—n), one obtains the remarkably simple result, 



One obtains therefore, 

0(i;x) = i J dnR{nf 



(3.11) 



independent of the the observation point x: the potential 
interior to a hi-homoeoid is constant. In a more New- 
tonian style, Chandrasekhar [3] demonstrates this geo- 
metrically by showing that the electrostatic force due to 
infinitesimal elements of charge at ±n precisely cancel 
one another. Using (3.4) and (3.5), the potential in the 
interior of any s-homoeoid is also constant and given by 

(j>{-K) = / d/xp(/i) = (/)(1,x)['(/;(m2) - ■0(Mi)], 

J 111 

(3.12) 

in which, for later convenience, the integrated charge den- 
sity function. 



V'(/i) - / p(/i')V 

JQ 



(3.13) 



has been defined, so that p = dip/dii. For a fi- 
nite homogeneous homoeoid, bounded by /i = /ig and 



/i = 1 and with uniform charge density po, one obtains 
ilj{p.2) - ipipl) = (1 - A*o)po in (3,12). 



C. Representation of the basic integral in terms of 
elliptic functions 

One may reduce the angular integral (3.11) to a more 
convenient form as follows. First, choose the usual spher- 
ical angular coordinates 9 and (j) to be centered at the ori- 
gin. In order to adapt the present notation to that which 
is standard for elliptic integrals, one uses oi as the z-axis, 
02 as the y-axis, and 03 as the x-axis. Chandrasekhar 
[3] makes a different choice, but the final answers must 
clearly be fully symmetric in the components of a. One 
obtains 



1 _ COS2(0) , ^. 2 



R{ny 



sin-' (6*) 



COS^((/)) 



sin (0) 



(3.14) 



Substituting this form into (f>{l;x.), and using dfJ = 
sm{6)d0d(f), one first performs the (j) integration [using 
the substitution u ~ tan((/))] to obtain 



(l;x) = 27ra^a2a3 



n/2 



sec^{9) sin{9)de 



yja\ + a\ tan2(6') v'ag + a\ toxoid) 



(3.15) 



r 



The substitution t — a\ tan^ (S) now simplifies this to 

</.(l;x) = i.(a)A(a), (3.16) 



in which 



A(a) ^ 
A(a,t) = 



dt 



A(a,t) 
^{al+t){al + t){al + t). (3.17) 



Prom Ref. [12], p. 220, with the convention oi > 02 > 
03, A{a) may be expressed as: 



A(a) 



=F((p,fc) 



(3.18) 



in which the elliptic integral F{lp^ k) is defined by 



^l-fc2sin'(a) 
dx 





,/{l-x^){l-k^x^)'- 



with parameters 



sin'^(iy9) = 1 — Og/o^ 



Oi — OS 



(3.19) 



(3.20) 



The result (3.18) follows directly from the substitution 
X = y/ (of — a|)/(o^ + 1) in (•J li ), or equivalently the 
substitution x — -^/l — ai/af cos(6') in (''- l^>). For later 
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purposes, we define also the auxiliary elliptic function 



D. Potential exterior to an hi-homoeoid 



dt 



A{a,t) al +t 



(3.21) 



The two functions obey 

aj ' dip 
E-F 
k 



dE 
dip 
dE 
'dk 



al dF _ al 



dF 



1 

k^ 



E-k'^F as al 
k a2\al 



- 1 



(3.22) 



Thus if one has available numerical algorithms to eval- 
uate E and F, then all derivatives of E and F follow 
by algebraic manipulations alone. This will be used be- 
low to generate an iterative procedure for computing all 
required integrals. 



The remarkable cancelation that produces { • . ) fails 
when the observation point x is external to the hi- 
homoeoid, which will be denoted Ei. To make progress 
one must use a different approach. The key idea is to seek 
the equipotential surfaces of (/)(l,x). One may reason- 
ably guess that such surfaces form a family of ellipsoidal 
surfaces, but it is not obvious what the corresponding 
family of principal axes should be. For a given x, let 
a' = (a']^, 02,03) label a concentric hi-homoeoid E2 pass- 
ing through x: X^a^a/'^a = 1- One seeks conditions on 
a' such that this ellipsoidal surface is an equipotential. 

For each point xi on Ei , let a corresponding point X2 
on E2 be defined by the relation X2^al o!^ = xi^a/a-a- This 
correspondence may now be used to map integrals over 
the surface Ei into integrals over the surface i?2. Volume 
elements then translate as d^X2 — {a'la^a'^l aia2a'i)d^xi. 
The potential at x may now be manipulated as follows: 



<^{Eu^)-l ^^fl-E%l=^/ "'"^ -^fl-E%l (3-23) 

where x' is the point on Ei corresponding to x: x'^/ua = Xa/a'^. If it were not for the altered Coulomb factor in the 
denominator, the second line would correspond precisely to the potential at a point on Ei due to a charged homoeoid 
E2- The idea now is to choose a' in such a way as to restore the Coulomb factor to the form |x' — X2I. One therefore 
seeks a solution to the equation, 

I 



If — = A is independent of a, the ellipsoidal condi- 
tions J2a^a/'^'i = 1 = J2a^2,a/C''a make the final sum 
vanish identically for any choice of the two points x' and 
X2 (or, equivalently, x and Xi). Two ellipsoids related by 
this condition are called confocal, and the corresponding 
"covariance" of the Coulomb factor is known as Ivory's 
Lemma [3]. 

One obtains therefore the following remarkable result: 
if A > labels the unique ellipsoidal surface E2 confocal 
to El and passing through the point x, i.e., 

y = 1 (3.25) 

then 

MEi,^) = ^ 0(^2, x). (3.26) 

^ ^ V(«^ + A)(a^ + A)(a^ + A)'^^ ^ ^ ^ 



The prefactor on the right hand side may be stated geo- 
metrically as a condition that Ei and E2 carry the same 
total charge. To conclude the argument, note that since 
(/)(i?2, x') is the potential in the interior of a homogeneous 
homoeoid, which by the previous results is a constant in- 
dependent of x', the original potential 0(i?i,x) will be 
a constant, independent of the point x on the confocal 
ellipsoidal surface E2 ■ This verifies that the family of con- 
focal ellipsoids, which sweep out all of the space external 
to the hi-homoeoid as A varies over the range < A < cxj, 
characterize completely the equipotential surfaces. The 
completeness of this family of surfaces also demonstrates 
that there can be no further independent solutions to 
(3.24). Furthermore, (3. 20), together with (3.17), yields 
the explicit formula 
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</.(£;i,x) =t;(a) / 
Jo 



dt 



y/{a\ + A + t)(a2 + A + t){al + \ + t) 



w(a) 



dt 
A 



u(a)A(a, A). 



(3.27) 



From (3.17) and (3.18) one then obtains the general el- 
hptic integral representation, valid both exterior and in- 
terior to the hi-homoeoid: 



0(^1 ;x) 



w(a) 



(3.28) 



in which k remains as defined in (.)._ ,), while (p(A) is now 
obtained from 



sin2[^(A)] = l-:^ 



at + A 



(3.29) 



E. 



Potential in the interior and exterior of an 
s-ellipsoid 



One may finally use the results (3.17) and (3.27) for the 
potentials inside and outside a hi-homoeoid to construct 
the full potential due to any s-ellipsoid, E. From (3.4) 
and (3,-j) one obtains for the exterior potential 



pl /"OO 

(x) = u(a) / d/ip(^) / 



dt 
A' 



X e 



(3.30) 



in which E"^ is the complement of E^ and A(/i) is the so- 
lution to the equation 2^a/["a + ^ A*: ^-i^d there- 
fore parameterizes the confocal ellipsoidal surface passing 
through the point x/y^. This formula includes the case 
of the general s-homoeoid which simply corresponds to 
a vanishing p(/i) for /i smaller than some /xq > 0. Inter- 
changing the order of the integrations, and noting that 
A(/i) — ?> oo as /i — > 0, one has 



(x) = w(a) 
= w(a) 



dt 



pi^J.)dfj, 



dt 
A 



mi) - l|;[^im xei?^(3.31) 



in which, as in (3.27), A = A(l) labels the confocal ellip- 
soidal surface passing through x itself, tpifi) was defined 
in (3.13), and fi{t) = E„ xl/{al + t). 

To calculate the interior potential, one divides the con- 
tributions into two parts: the potential from the 
similar s-ellipsoid whose surface passes through x, and 
'/'out from the complement s-homoeoid whose inner sur- 
face passes through x. From (3.31) one obtains 



i(x) = w(a) / dfip{fi) 

Jo "'A(ai) 

= «(a) r^{^[^(0)]-^[MO]}, (3.32) 



°° dt 
A 



where /i(0) = J^a ^a/'^a labels the surface of the interior 
s-ellipsoid. The lower limit on the t integral vanishes 
because x lies right on this surface. On the other hand, 
from ( ) one obtains 



0out(x) = w(a) 



^ f°° dt 

°° dt 
A' 



/m(o) ^0 
{^(1)-^[M0)]} 



(3.33) 



The total interior potential is then given by 
dt 



0(x) = v{a) 



{V(l)-^[/i(i)]}, xei?. (3.34) 







The only difference between (3.-)] ) and (3.34) is the lower 
limit A on the t integral. 

As a simple application of these formulas, for a uni- 
formly charged ellipsoid, p{fi) = po, one has V(l) ~ 
■)Jj[pit)] = po[l - J2a ^l/io-l + and therefore 



(x) = pov{a) 



A(A)-5]4aW(A) 



in which 



4^' (A) 



dt 



For interior points one simply sets A = 0. 



(3.35) 



(3.36) 



F. Generalization to certain classes of 
non-s-ellipsoids 

The monomial charge density p(x) = x^^ . . . x^^ corre- 
sponding to (2.45) clearly does not fall into the s-ellipsoid 
category. The following trick, however, may be used to 
generalize the results to include this form. Given any 
solution to the Poisson equation — V'^(/)(x) = 47rp(x), 
by taking derivatives of both sides it is a trivial ob- 
servation that da4>{x) is the solution for charge density 
9qp(x). Charge densities of the form p = p{p) with 
p = ^"al^^ have so far been considered. One obtains 
therefore 



dpp 



dpd^p 



"■fi 

III x4a;^a;-j, 
^ ^^^^^ 



(3.37) 



-Si 



and so on. Repeated applications of the derivative trick 
therefore produce s-ellipsoid charge densities multiplied 
by polynomials in the components of x. Clearly (2.45) is 
a special case of this. 
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IV. COMPUTATION OF THE BASIC 
INTEGRALS 

All the necessary tools to compute the potentials 
2^m(x), equation (3,1), are now at hand. Motivated by 
(3.37), before specializing to these, we consider the more 
general class of integrals 



P(^)(x)= / d 
Jv. 



( f\ f ''^1 f ''^2 / 



(4.1) 



in which q = q{fj.) only. These reduce to (2.45) when 
q = 1. The formalism developed in this section allows 
analytic treatment of all integrals of the form (4.1) by 
selecting of an appropriate p(/i) upon which to apply the 
derivative trick. The result for Vq (x) is contained already 
in (3.31) and (3.34): 



2?o(x) 



via) 



|X — X' 

°° dt 



pMt)] 



(4.2) 



in which A(x) = for x interior to the ellipsoid, and is 
determined by (3.25) for x exterior to the ellipsoid, and 



(4.3) 



Thus q = —dpi/dfi and both q and piifi) vanish for fi > 
1. For future reference, the sequence of charge densities 



Pnip) are defined iteratively via 

Pn+l{p) = Pn{p')dp.', 71=1,2,3, 



(4.4) 



so that pn = —dpn+i/ dp,^ and all pn vanish for /i > 1. 

Next, from the first line of (3,37) one has qxa 
— {a'^/2)daPi, and hence 



= -~\alda J d^x 



3 / PUX 



f°° dt 



in which ei — (1, 0, 0), etc., and a term coming from the 
x-dependence of the lower limit of integration A vanishes 
since ^(A) = 1, and hence p2[p{\)\ — P2(l) = 0. 
A similar procedure, using the identity 



qxaxp = ^a^aldadi3P2 + ]^a\pi5ai3, (4.6) 



yields 



2?e„+e3(x) = j d-'x 



-al^alv{a.)xax,3dadfs / -^Pz[p{t)] + -a^w(a)(5a/3 / -irP2[p{t)] 



[°° dt 1 [°° tdt 

= v{a)alalxaxp XT-^— TV^^— 7YPi[M(i)] + o"(^)"a'^"/3 / 2 ..-. P^ipi*)]- (4-7) 

Jx ^iai+t)[ap+t) 2 A{ai+t) 

Similarly, the identity 



VaXisx^q = - ^alaja'^dadpdjps - j [a\al^5af3d^P2 + (7 O /3) + (7 O a)] , (4.8) 



yields 



•D6„+6^+e^(x) = / d-^x 



,3 ,l{^')^'a^'l3^'j 



X — X' 



f°° dt 
v{.)alay^x^x,x, A{al + t^a} + t){a^^ + tf'^^^*^^ 

1 f /"oo tdt 

+ 2^(a) ^ala^^x^S^p A(a2 + ^)(a2 + ^) ^'^^^^^^ + ^ P) + {1 ^ a) \ . (4.9) 



In all cases, terms arising from derivatives acting on the lower limit A do not contribute since Pn{l) = for n > 1. 
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It is clear that calculations become substantially more 
tedious as n increases. General cases beyond n = 3 will 
not be treated explicitly in the present work. Simplifica- 
tions for solid homogeneous ellipsoids that allow general 
evaluation of I?m(x) will be treated in Sec. V . 



It is apparent that interior to the ellipsoid, where A = 0, 
(f)n is an even polynomial of degree 2(n + 1). 



A. Monomial charge densities: even cases 



V. EXPLICIT EVALUATIONS FOR THE CASE 
OF A SOLID HOMOGENEOUS ELLIPSOID 

For the special case that the scatterer is a solid, ho- 
mogeneous ellipsoid, q(x) = 1, which will be the focus of 
all explicit computations in this and later sections, one 
obtains Pn{fi) = (1 — for < ^ < 1, p„(^) = for 

/i > 1. From ( ), the potential due to p„ is given by 



0n(x) 



1 

nl 



V, |x-x' 



dt 



(« + l)!A A 



1-E 



al+t 



n+l 



(5.1) 



-'(a)i: ,;^L;:.;;^^i,/ .(a.A). 

in which we define |k| — and the prime on the 

last sum indicates that it is restricted by |k| < n + l, and 
we have defined the basic integrals 



.(a,A)= £ 



dt 



A(a,i)nj«^+t)* 



(5.2) 



with Aq = A [see (3.17)]. These obey the simple iterative 
relation: 



d 



-4k = -(fc„ + l/2)Ak+e„, 



(5.3) 



in which all derivatives are at fixed A. Thus, 

(-i)i''ir(i/2)3 



Ak(a,A) = 



r(A:i + l/2)r(fc2 + l/2)r(A:3 + 1/2) 

aiki 

^A(a,A). (5.4) 



9(a2)fcia(a2)fe2a(a2)fe3- 



One may now use (^ ) to derive formulas for the po- 
tential due to monomial charge densities, extending the 
results of Sec. 1 \ to general index values. From the first 
line of (5.1), one obtains for even monomials 

I52i(x)^/ ^!^[](xL)2'. =ai0|,|(x;a) (5.5) 
Jy, |x-x'| 

where, again, |1| = la, and the operator acting on the 
a-dependence is defined by: 



5i^(-l)l'l 



am 



(5.6) 



The derivatives eliminate all monomials of order lower 
than |1|, and there are no contributions from the a- 
dcpcndcncc of the integration region Vs because the inte- 
grand of (pn, along with its first n — 1 derivatives, vanish 
on the boundary. This observation allows us to apply the 
same derivatives to the last line of ( /. i) at constant A: 
terms arising from derivatives of the a-dependence of A 
must, via the second line of (" ' cancel. One obtains, 
therefore, 



2k2 2fe3 
''2 ^3 ' 



(5.7) 



in which the prime on the sum indicates the constraint 
|k| < |1| -I- 1, and the coefficients are given by 



(-1) 



|k| 



fci!fc2!fc3!(|l|-|k| + l)! 



5i[^;(a)Ak(a,A)]. 



(5.8) 

Internal to the scatterer, v'^\x) is a polynomial of de- 
gree 2(|1|-H1). 
By iterating the relation 



one obtains 



-,2m+l 



Y = fc + m + - 



(-1) 



d a 



2m+l 



-^2m+5 



2/ (a2+t)fe+™+5 



— m + - 



-,2m+3 



2/ (a2 +t)fe+™+5 



2(m+;+p)-|-l 



- = (-]]'^k^ niklm) _ 



in which we have defined the coefficients 



k\ \P J T{\~-ra-l) V{k + m+\) 



(5.9) 



(5.10) 



(5.11) 



14 



This produces the expUcit form 

= (111 -t\+ i), E'^'-+p(^.^)ng£°'"°^"^^^°^ (5-12) 

villi ' p a 

where the prime indicates that the sum is hmited to the range < < Inserted in to ( ), the result (•") I'i) 
exphcitly exhibits the monomial integral {^>-'>) to a linear combination of the basic integrals (■^•-*)- 

If one wishes to compute the coefficients V^^^ for all (even) monomials of degree at most A'', since < |1| < [N/2], 
and hence < |k| < [N/2] + 1, one therefore needs to compute all basic integrals with indices constrained by 
< |n| < 2[N/2] + 1. 



B. Monomial charge densities: odd cases 



Consider now cases where the monomial is not even. The desired results follow immediately from the identities 



5i92C'30n+3(x) = 





f d^x' 




|x-x'| 




f d^x' 


n\i 


|x-x'| 




f d^x' 




|x-x'| 



1-E 

7 



a ^ 



(5.13) 



By applying the operator d\ to these expressions one may 
now isolate the required monomials: 



^21+.n(x) = (~l)l'+™lai 



(5.14) 

where the elements of m are all either or 1, and as 
before |1 + m| = ^Yli-^^a + Wq.). Following the derivation 
of ("(. l ^), one therefore obtains the monomial expansions 



2k„ 



(5.15) 



in which the sum continues to be restricted to |k| < |1| + 1 
and 



D 



(21+m 

2k+m 



'(A) 



Ill-Ikh 



"YyjE ^k+m+p(A) 



(ia+Pc+m„) 



(5.16) 



where the sum is again over all < Pa l£ la- Equation 
(5.12) is now clearly a special case of (5.16) in which m 
vanishes. If one desires these coefficients for all monomi- 
als of degree at most N, then the basic integrals A„ will 
be required for all n such that < |n| < + 1. 

It is straightforward to check that (5.15) and (5. 10) 
reproduce the partial results in Sec. IV if one sets q = 1- 



C. Recursion relations for the A integrals 

We finally reduce the computation of the A integrals 
to algebraic recursion relations, given only the pair 



^000 — 



A,oo = - ^ (5.17) 

[see (3.18)-(3.22)]. 

The relations are based on the following three identi- 
ties. First, by integrating the identity 



d, 



1 1 sr^ h + ^ 



over the range A < i < 00, one obtains the relation 

1 



E + J)^k+e,(a,A) 



(5.18) 
I 

r. (5.19) 



Second, it is easy to verify the homogeneity relation 
Ak(Ka, K^A) = K-2(fci+fc2+fc3)-i^j^(a^ x). (5.20) 

for arbitrary scale factor k > 0. By taking the derivative 
of both sides with respect to k and setting k = 1, one 
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obtains the Euler relation 



{al + X) (^fc^ + ^k+e^ (a, A) 



|k| + - ) Ak(a,A). 



(5.21) 

By combining (^i Hi) and (' ' ' ' ) one obtains for each 7, 
(«7 - ) f'^" + 1) ^k+e„ (a, A) (5.22) 



a^ + A 



|k| + i)Ak(a,A). 



nj«^+A)^-+^ 

Third, by noting the trivial identity 

(al + t)- (aj + t) (5.23) 
one obtains for any a 7^ /3, 

Ak(a, A) = -^-^ — 2- [^k-e„ (a, A) - ^k-e^ (a, A)] . 

(5.24) 

Equation ( ) provides a relation between and 

any two of its "forward" neighbors, while (5.24) provides 
a relation between ^k and any two of its "backward" 
neighbors. It is easy to see that by applying these two 
relations appropriately one may recursively generate any 
Ak from the pair (5.17) alone. In numerical implementa- 
tions care must be taken, however. The recursion (5.24) 
will likely be unstable if \aa — <if)\ is too small. In this 
case an approximate form for ^k will need to be com- 
puted for large enough |k|, and the recursion relations 
iterated backward to smaller values [14]. 

VI. CASES OF DEGENERACY 

Great simplifications occur when two or more of the 
semi-major axes are identical. The cases of most interest 
are spheroids where ai — a2 = a, and either > a 
(prolate spheroid) or 03 < a (oblate spheroid) . The third 
case, oi = 02 = as = a, corresponds to the trivial case of 
the sphere. 



A. Sphere case 

For the case of the sphere one finds the fully analytic 
result 



dt 



X (a2+<)^+i 
1 1 



(6.1) 



where K = |k|. From (3.25), for exterior points one 
obtains 

A = |x|2 - a^, |x| > a. (6.2) 
B. Prolate and oblate spheroids 

For the more interesting case of spheroids, the basic 
integral 



A(a,A) = 



dt 



is elementary, and for the prolate case one obtains 



(6.3) 



A(a, A) = 



2 _ „2 



:ln 



ai + \ 



ai + \ 



(6.4) 



while for the oblate case one obtains 



^(a, A) 



\/a2 - aj 



arctan - 



a? + A 



(6.5) 



In either case, the higher order integrals 

AkkJa,X)= 7- — ., „ — ,1 , (6.6) 
Jx (a2 + t)'^+i(a^ + <)'=3+2 

where k — ki + k2, may now be generated by iterating 
(5.22), which may now be put in the form 



Afe+i,fe3(a, A) 
^fe,fe3+i(a. A) 



1 



1 



1 

{ks + ^)ia^-al) 



(a2 + A)fc+i(a2 + A)'=3-i 
1 

(a2-f A)fc(a2-HA)^3+i " 



^ + ^3 + 2 ) ^kks (a. A) 



k + k3 + - )^fcfc3(a,A) 



■7) 



Here, the first and second lines follow by setting 7 = 3 and 7 = 1 or 2, respectively. Note that care should be taken 
for small jaa — a\ since the terms in brackets are then small as well, so that the overall result remains finite. 



16 



Finally, equation (3.25) for A is quadratic, with solu- 
tion 



- -(a^ + a^-|x|^), 
which is consistently positive for x outside Vg- 



VII. FORMAL SOLUTION TO THE 
SCATTERING PROBLEM 

Now that a complete formalism for the evaluation of 
the necessary integrals has been presented, one may turn 
finally to the solution of the scattering problem (2.'l) us- 
ing the high contrast formulation (2.24). The main work 
in the application of the mean field approach is the com- 
putation of the electric field internal to the scatterer, 
which involves the evaluation of the I?™(x) [equation 
(2.4G)] inside the ellipsoid. For this purpose, one may 
therefore set A = in all of the formulas derived in Sees. 
V and (VI). 



Thus, the background field E(, 
noninductive, and one identifies 



-V$f, must also be 



$6(x) 



,n(r') • V'$(r') 



47r|x - 



47r|x - 



(7.3) 



Note that a solution to the Laplace equation is uniquely 
specified by a Neumann boundary condition, so there is 
no contradiction in the fact that the full $ appears in the 
second line, whereas only its boundary normal derivative 
appears in the first. 

If $ is a polynomial, the last integral is precisely of the 
type considered in Sec. . The two gradients imply that 
is a polynomial of the same order. In particular, the 
linear form $6(x) = Xa leads to 



$(x) 



Xa 



2n 



k2 v(a)Ae„(a, 0) 



(7.4) 



This is equivalent the well known solution for an ellipsoid 
in a uniform static applied electric field Eq: 



Eint — E] ^aXaEo^a- 



(7.5) 



A. Noninductive solutions 

Before proceeding further, it is worth further clarifying 
the nature of the projection in ( r , ), which removes the 
noninductive part of the electric field. This is most eas- 
ily done directly from the integral equation (2.0), where, 
for simplicity, we consider the case of homogeneous back- 
ground as well as scatterer. Using translation invariance 
of g = .9(|x ^ x'l) [equation ( )], and integrating by 
parts, as in (^.i-i), one obtains 

E(x) = Efc(x) + / d^x'g{x,x')E{yi') 
Jv^ 

/ dVg(x,r')n(r')-E(r'), (7.1) 



where Q 



nr and we have used V • E = 0. The 



second line is a pure gradient, and hence part of V^. In 
(2.22) we restricted the class of solutions to those with 
vanishing normal component. However, it is clear from 
this term that a contribution with n • E = 0{k^) [which 
provides the leading finite Kb correction to the boundary 
condition (2.22)], although contributing negligibly to the 
inductive part of E — which is contained entirely the first 
line of (~ J ) — does make a finite contribution to $. 

To explore this further, let us seek purely noninductive 
solutions E = — V<f>, which requires also = 0. This 
separation is consistent only to leading order in k^/k^, 
where (7.1) reduces to 



Efa(x) = - 



,n(rO ■ V$(rO 
47r|x - r'l 



(7.2) 



Higher order polynomials may be handled in a similar 
fashion. For example, the form <!>;, = XaXp (a ^ /3), leads 
to 



$(x) 

Xq/3 



^2 



k2 u(a)(a2 -f a2)^g^+g^(a, 0) ' 



(7.6) 



A somewhat messier calculation allows one to compute 



the response to 



Co 



2 



(a ^ /3) in the form 
with coefficients obeying the con- 



straint — 0. Progressively higher order polynomial 

forms for may be reduced to progressively higher or- 
der matrix inversions for the polynomial coefficients of 
[15]. 

It is clear from ( ) that all such solutions generate 
very smah internal fields, |Eiiit|/|Eb| = 0{h1/k'^). How- 
ever, via the last term on the right hand side of ( ) , 
the external field is of the same order as the background 
field, but, once again, is invisible to a purely inductive 
measurement. 

One may also compute corrections to these solutions, 
in powers of k^/k^, arising from the terms in (7.1) left out 
of {~i ■'■>)■ For metallic targets, however, these corrections 
are negligible, and one concludes that the noninductive 
modes respond quasistatically in the frequency regime of 
interest here. The physical origin of this result is that 
these modes are controlled by the induced polarization 
charges on the surface of the target, which respond essen- 
tially instantaneously to the applied field. Current flow, 
required to generate magnetic fields that would influence 
an inductive measurement, is negligible. 
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B. Perturbation theory for inductive excitations 

Consider now solutions to the internal field equa- 
tion (2.24). Simple dimensional analysis shows that the 
Coulomb integral is of order 



kjjLcrLil 



2t:51 



(7.7) 



where Ls is target characteristic diameter, and Ss{uj) ~ 
c/^/2TTcrfiuj is the target skin depth [1]. The parameter 



Vs = 



Ls _ nLs 
(5. 5 cm 



1/2 



W S/m 



1/2 



/ 



1/2 



100 Hz 

(7.8) 

[written here in MKS units; compare (2,iV)] is small at 
low frequencies, and is the formal mean field expansion 
parameter. The dimensional quantities inserted here typ- 
ical of compact metallic targets that might be of inter- 
est. To zeroth order, for ry^ ^ 1, equation (2.1!)) reduces 
simply to A = Af,: the target is transparent to mag- 
netic excitations [as opposed to its essentially complete 
opaqueness (7.3) to polarizing fields]. In particular, the 
zeroth order solutions 



''Imp 



(x), 



(7.9) 



may be parameterized by the basis functions defined by 
(2.39), and rescaled according to (2.32). Here, the V$6 
part of Ef,, which has been projected out of (7.9), gen- 
erates, via Sec. VHA, 0{kI/k^) polarization corrections 

to Eint. 

One may conveniently compute corrections iteratively: 



E(x) = f:f^)"AE„(x), 
in which the terms satisfy the recursion relation 



(7.10) 



Magnetic field correction AH^ : a^/a = 
H along z 



0.5; X = plane 





FIG. 1: (Color online) Leading magnetic field correction 
(T.i")) in the yz-plane for a spheroid, ai = a2 = a, using as- 
pect ratio as/a = 0.5. Upper plot: background field (~ ' ') 
taken along the z-axis. Lower plot: background field taken 
along the j/-axis. In both cases, for a; = the field also lies in 
the yz-plane. For background field along x (not shown), the 
field in this plane is along x as well. 



AE„+i(x) ^T'. / d'x 



,3 ,AE„(x) 



(7.11) 



beginning with AEo = Va^b- The series (7.10) is 
an explicit expansion in powers of the small parameter 
kfjLa/c = 77g/27rLj. Since the basis functions (7.9) are 
polynomials of degree N = l + 2p + i — l, and application 
of the Coulomb integral adds two to the degree of the 
numerator (application of the projection operator does 
not change this), the nth term in the series (7.10) is of 
degree N -\- 2n. 

Consider, for example, the linear form 



AEo = —Xaen 



-xpea, {a ^ 13), 



(7.12) 



which circulates in the a/3-plane, and indeed lies in the 
space of functions defined by (2.22). These are, in fact, 
three independent linear combinations of the three / = 1, 



P 



harmonic basis functions z|j^p^ 



This form corre- 



sponds to a constant magnetic field 



V X Eb 



1 < 



ikyL aaCLp 



(7.13) 



where (ck/37) is taken to be a cyclic permutation of (123), 
and the magnetic field is therefore orthogonal to the a/3- 
plane. This will be an adequate model of the background 
field if the transmitter is sufficiently far from the target. 

Applying the results of Sec. \ [specifically, ( ) and 
the first line of (' ') with n = 1], one obtains the leading 
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correction 
AEi — aQa^w(a) 



(a, 0) - ^ Ae„+e^ (a, 0)xl 



(7.14) 



The projection operator (especially applied to the cubic 
terms) is very messy, and its result will not be displayed 
here. However, since it subtracts only a gradient, the 
curl of (~ ' 4) directly produces the leading correction to 
the magnetic field: 



AHi 



V X AEi 



gga/jf (a) 
ikfi 



(eT,{^e„(a,0) +^g^(a, 0) 

,(a,0)+3A2gJa,0)]a:2 
[^e„+e^(a,0) +3yl2§^(a,0)]x^ 

[Ae„+e.,(a, 0) + Ae^+6^(a,0)]4} 

2x^ [Ae^+e-,{a,0)xaea 



A, 



,(a, 0)a;^e^]), 



(7.15) 



in which, again, {aj3^) is cyclic permutation of (123). It 
is easily verified that this correction is divergence-free, 
as required. Example field patterns are shown in Fig. ' 
for the case of an oblate spheroid (discus-shape), where 
the required A-coefficients are computed in closed form 
using the results in Sec. ' " . 

It is clear that the complexity of the polynomial form 
AE„ increases rapidly with n, and a numerical imple- 
mentation is required to keep track of all the terms. In 
later sections we will show theoretical results, and com- 
parisons to experiment, using all 232 basis functions with 
I + 2p < 7, where convergence is found even when rjs is 
not small. 



C. External field 

It is the field external to the scatterer that is relevant to 
target detection. Equation ( .: ) allows computation of 
(the inductive part of) the external field once the internal 



field is known [16]. Specifically, from the series (~ IH), 
one obtains 



e: 



('"'*)(x) = zfcA(x) = ^ 



ikfia 



AE(,"'i)(x) (7.16) 



in which AEq"'^'' = ikAh is the inductive part of the 
background field (both inside and outside the target), 
and 



AE(^'"'^)(x) 



d^x'^^^^. 



K, n>l. 



|x - X' I 

(7.17) 

If the result for AE„_i is a polynomial, as described in 
Sec. VII B, then the results of Sec. V directly produce the 
external field in the form of sum of products of polynomi- 
als with elliptic functions. The latter are now functions 
of position x via the nonzero value of A(x) — see (3.25). 

For example, the linear form ( ) for the background 
field produces a leading correction 



^j,(ind) ^ aaapv{a) 

X < Xaep 



xpt 



Ae^ (a. A) - ^ ^e„+e„ (a, X)xl 
(a. A) - ^ ^e„+e^ (a, A)a: 



(7.18) 



which is identical in form to (7.14), but with the projec- 
tion operator omitted, and now including nonzero A. 



D. Far-field asymptotics: multipole expansion 

Computation of the external field greatly simplifies if 
one is far from the target — greater than a few maximum 
target radii, say — but not so far that non-quasistatic cor- 
rections to the Green function (2.4) become important. 

First, for large |x|/Ls the cubic equation (3.25) for A 
has the expansion 

A(x) = |x|2_^a^x2^0(l/|xp), (7.19) 

a. 

in which x = x/|x| is the unit vector with components 
Xa — Xal\y^- For large A one finds in turn. 



Ak(a,A) = 



dt 



dt 



t|k|+3/2 



l-^E(2fca + l)ai + 0(1^2) 



Ej2/c„ + l)al 1 



2|k| -h 1 Alk|+i/2 



2|k| 



Alk|+3/2 



o 



2t 
1 



Alk|+5/2 



(7.20) 




FIG. 2: (Color online) Spheroid decay rate spectrum vs. aspect ratio a = a-i/a, with conductivity a and radius a fixed as 
indicated. Unit permeability, ^ = ^f, = 1, is used as well. The first 125 decay rates (out of a total of 232 that are computed 
using all 232 basis functions with order N = I -\- 2p < 1) are plotted for each 0.1 < a < 10. Apparent sharp bends in the 
curves at the top of the spectrum are artifacts of this truncation. Blue dots at a = 1 show exact analytic results for the sphere, 
with degeneracy effects from enhanced symmetry evident. Red dots a.t a — 4 mark the four (doubly degenerate) vertically 
circulating modes shown in Fig. blue circles mark the four azimuthally circulating modes shown in Fig. I. The green dot at 
the right is the analytic result for the lowest mode for an infinite cylinder, a — >■ oo. More slowly increasing branches at the left, 
Q! ^ 1, correspond to modes with current patterns circulating in the xy-plane. The red dots at q = 0.1 mark the first twelve 
such modes shown in Fig. 5. 



with, as before, |k| = J^a ^a- Substituting (^ H)), one obtains explicitly 



Ak{sL, A) 



1 



2|k| + l|x|2|k| + l 



E2 -2 



2|k| +3 



1 



|x|2|k|+3 



o 



1 



|x|2|k|+5 



(7.21) 



As an example, leading behavior of the external field 
correction (7.18) takes the simple magnetic dipole form 

(ind) 4,aaapv{a.) , _ 
AEi = ^^i^p ixaBf! - xpe^) . (7.22) 

Higher order terms in ( i iO) will generally all have a lead- 
ing l/|x|2 term, but with more complicated angular de- 
pendence. This may be formalized via a vector multipole 
expansion using the vector harmonics (2.35) [1]. 



VIII. TIME DOMAIN RESPONSE 

Having presented the general theory, and some leading 
order perturbation results in the previous sections, we 
now turn to more sophisticated applications of the the- 
ory that require numerical implementation of high order 
expansions in r]s- Specifically, theoretical predictions for 
time-domain EM measurements will now presented (see 
Sees. IB and II D) and compared to experimental data. 
As described previously, time domain measurements are 
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FIG. 3: (Color online) Four lowest order vertically circulating mode shapes, and corresponding decay rates (red dots at aa/a = 4 
in Fig. ), for a 10 x 10 x 40 cm radius aluminum prolate spheroid: {Ey, Ez) are plotted in the x = plane (where Ex vanishes 
identically) . 



Azimuthal modes: E (x=0) 

<l> 

Mode 3; ?i= 15.7 s"'' Mode 8; A, = 20.6 s"' Mode 9; A, = 25.6 s"^ Mode 14; X = 33.7 s"' 




-10 10 -10 10 -10 10 -10 10 

x-coordinate (cm) x-coordinate (cm) x-coordinate (cm) x-coordinate (cm) 



FIG. 4: (Color online) Four lowest order azimuthally circulating mode shapes, and corresponding decay rates (blue circles at 
as/a = 4 in Fig. 2), for a 10 x 10 x 40 cm radius aluminum prolate spheroid. The azimuthal component of the electric field 
E^ (or Ex), normalized by its maximum magnitude, is plotted in the x = plane; all other components vanish. Currents in 
the lowest mode (far left) circulate in the same direction at all heights. As the mode order increases (from left to right), the 
number of nodal surfaces, across which the circulation direction of the current changes sign, increases. 
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Oblate spheroid drumhead modes: 2D plots of (E ,E ) at z=0 cm; aspect ratio = 0.1 
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FIG. 5: (Color online) Twelve lowest order horizontal circulating "drumhead" mode shapes, and corresponding decay rates (red 
dots at as/a = 0.1 in Fig. ), for a 20 x 20 x 2 cm radius aluminum prolate spheroid: {E^, Ey) are plotted in the z = Q plane 
(where E^ vanishes). Modes 1,6,15 have azimuthal symmetry (angular momentum index m = Q) and are non-degenerate. The 
remainder, with m > 0, are doubly degenerate, with the second mode shape obtained via a ■n/2m rotation about the z-axis 
(modes 2,9 have m = 1; modes 4,13 have m = 2; modes 7,18 have m = 3; mode 11 has m = 4; mode 16 has m = 5; mode 
23 has m = 6). For fixed m, the higher order modes have increasingly complex radial structure, with the circulation direction 
changing sign with radius (see modes 6,9,13,15,18). 



those of the freely decaying response of a target after 
termination of an applied pulse. 

Any time-domain response may, of course, be writ- 
ten as the Fourier superposition of a spectrum of fre- 
quency domain responses, each of which may be individ- 
ually computed using the previous theory. However, as 
will be described below, the formulation in terms of a 
superposition of freely decaying modes, equation (11), 
is numerically more efficient because these need only be 
computed once for a given target, avoiding recomputa- 
tion of the perturbation series for each of a continuous 
set of frequencies. In fact, these modes can be used to 
directly solve the frequency domain problem as well. 

A more important point is that a rapidly terminated 



pulse (over tens of microseconds, in the experiments to 
be analyzed further below) has a spectrum that includes 
some very high frequencies (e.g., tens of kHz), for which 77 
is a far larger than can be handled at any achievable order 
in perturbation theory. However, this part of the spec- 
trum mainly excites very rapidly decaying modes that 
disappear from the later-time domain response. Thus, 
the mode approach allows accurate prediction of the sig- 
nal at later time even when it fails at earlier time. In fact, 
the very large number of modes appearing at early-time 
make it a poor representation of the response. A comple- 
mentary description in terms of the inward diffusion of 
screening surface currents [7, 8] should be used instead. 
It will be seen below that a merging of the two theories 
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provides a complete description of the signal. 

A. Mode computation 

The freely decaying mode computation is implemented 
via the generalized eigenvalue equation (2.4.3). Specif- 
ically, the modes are written as a superposition (2.42) 
of the truncated set of vector harmonic modes (2.39), 
rescaled according to ( ■), with I + 2p < N ioi some 
chosen upper limit N. Results will be shown for N = 7 
(which yields a total of 232 basis functions — 116 for each 
type i = 1,2), which will be seen to suffice for accurate 
comparison with experiment. Since zj^^ (x) is a poly- 
nomial of degree Z + 2p + i — 1, the x' integral defining 
the iJ-matrix (. ) produces, via the results of Sec. 
a polynomial of degree l-\-2p + i + l. The final x integral 
[which, it should be recalled, automatically implements 
the projection operator Va in (2,24) and (2.41)], in both 
(2.29) and the definition of the 0-matrix (2.28) (with 
uniform a here) is then a pure polynomial integral, and 
is trivial. 

Figure 2 show results for the decay rate spectrum of a 
range of spheroids (ai — a2 = a). The lowest 125 decay 
rates (out of a total of 232 computed at this at order 
N = 7) are plotted for each of 201 aspect ratios a = 
a^/a in the range 0.1 < a < 10. Even though these are 
plotted for a particular choices of radius and conductivity 
(as well as permeability), in the high contrast limit the 
combination A/icra^ is independent of all three. Hence 
these results may be trivially rescaled to obtain results 
for any spheroid with the same geometry. The mode 
eigenf unctions scale trivially with a as well. Note as well 
that the azimuthal symmetry (which is exactly preserved 
in by the perturbation theory at fixed N) means that 
the z angular momentum index m is a "good quantum 
number" , making the matrices O and H block diagonal. 
The eigenvales for ±to are also degenerate, so that many 
of the lines in Fig. - actually represent pairs of modes. 

Exact results for the sphere, at a = 1, are shown by 
cyan dots, and indicate the accuracy of the method for 
rather large effective values of rjs = 0(10) [obtained by 
substituting a for Lg and A for / in ( )]. The total 
angular momentum index I is now also a good quantum 
number, and the vast increase in degeneracy is evident. 

The infinite right circular cylinder of radius a corre- 
sponds to a — >■ oo. In this limit, translation invariance 
implies that the z-dependence of the modes is given by 
gzfez aj-i^ii^j-aj-y wavenumber k. The analytic result for 
one of the fc = modes (with current traveling up one 
side of the cylinder and down the other) , is shown as the 
green dot, and is seen to have converged even at a — 10. 

A few mode shapes for the prolate spheroid with a — 4 
are illustrated in Figs. 3 and 4. The four modes in Fig. 3, 
consisting of an increasing number of vortices circulating 
in the yz-plane and indicated by the red dots in Fig. 2, 
are doubly degenerate (with the second mode obtained 
by 90° rotation about the z-axis) . One may think of these 



modes as approximating those of an infinite cylinder with 
k w nir/a^, n = 1,2, 3, 4. 

The modes shown in Fig. , indicated by the blue cir- 
cles in Fig. 2, are non-degenerate. Here the current al- 
ways circulates in the xy-plane, but the flow direction os- 
cillates with z. Again, they may be thought of in terms of 
those of an infinite cylinder with k w wk j aj,, n — 1,2, 3, 4. 

A set of twelve modes for an oblate spheroid, with as- 
pect ratio a = 0.1 are shown in Fig. As seen Fig. 
2, since the radius a is fixed, the mode decay rates in- 
crease as they become vertically "squeezed" by decreas- 
ing 03. We denote these "drumhead modes" because 
a plot of the vertical component of the magnetic field 
Hz — dxEy — dyEz would look very similar to the surface 
height pattern of a vibrating drumhead. Of course, the 
latter oscillate at fixed frequency, rather than relax ex- 
ponentially, but there are parallels in the physics of the 
underlying mode patterns. 

There are two independent time scales that operate 
to determine the scaling of A with for a particular 
mode pattern, and the interplay between them basically 
determines the shapes of the curves in Fig. 2. 

For a circulating current vortex that is very thin (the 
dimension as for the modes in Fig. 5) compared to its 
horizontal extent (the dimension of one of the in- 
dividual vortices in the mode current pattern), the de- 
cay is dominated the decay of stored magnetic energy 
via Joule heating. Thus, the dissipated power scales as 
P ~ PR ^ (Jaya^)^ /{aa^) = J^a^a^/a, where J is the 
characteristic current density. The magnetic energy is es- 
timated as Uh — I^B'^ycSi in which B ^ I /cay — Ja^jc 
is the magnetic field within a circulating current field 
and VcS ~ is the effective volume over which it is sup- 
ported (mostly outside the target when a^/ay ^ 1), then 
scales as Uh ^ Pa^a\l(? [17]. The ratio A PjUn ^ 
(?la\iaa-i — / a ^o?a, as claimed. For the modes in 
Fig. 2, Otj/a remains roughly fixed as 03/0 — >■ 0, and this 
produces the observed A ~ 1/a scaling for a ^ 1. 

On the other hand, if two oppositely oriented currents 
streams lie very close to each other, the decay is dom- 
inated by transverse diffusion of the streams into each 
other, which causes them to cancel out. The diffusion 
time scale is tjj ^ DdP, where d is the current stream 
separation. 



is the diffusion constant, and one estimates A ~ 1/td ~ 
(? Iliadi}. For modes with alternating current sheets in 
the vertical dimension one has d oc 03 (e.g., the mode 
geometries pictured in Figs. '! or 4, where dj a-i ~ n for 
n — 1,2,3,4, but squeezed in the vertical dimension), 
and this leads to the A ^ 1 jct^ scaling seen for the steeper 
curves on the left side of Fig. 2. 
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Tx1 2: multiple bipolar pulse sequence 
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FIG. 6: (Color online) Transmitter pulse train generated by 
the NRL-TEMTADS platform. The pulse train (with maxi- 
mum current approximately 5.7 a), is bipolar, meaning that 
it consists of identical pulses of alternating sign. The average 
current is then zero, which helps average out certain tran- 
sients. The leading edge of each pulse is actually a sequence 
of three exponential relaxations, with time constants of 2.5 
fis, 0.33 ms, and 4 ms. The trailing edge is a linear ofFramp 
over a 10 /is interval. The total length of both the pulse and 
the following quiescent detection interval is 25 ms. 

B. Comparisons with experimental data 

1. Mode excitation and voltage amplitude computation 

Having explored some details of the mode physics, we 
now illustrate, through comparisons with experimental 
data on artificial aluminum spheroids, how the mode the- 
ory is used to compute measurable quantities through 
incorporation of a detailed model of the measurement 
platform. 

The first step is to relate the excitation coefficients 
An in the free decay series (11) to the source excitation 
pulse, encoded in S on the right hand side of the wave 
equation (2.3). In the time domain, this equation reads 

1 Attu 
— -atE(x, + V X V X E(x, t) = f 5Js(x, t). 

(8.2) 

One seeks the solution in the form 

E(x,t) = ^A.(i)e^"Hx) (8.3) 

n 

generalizing ( . > ). Substituting this into ( ), and using 
the mode defining relation (: ) and the orthogonality 
condition (2.4'i), one obtains the equation of motion for 
the individual amplitudes: 

(dt + \n)An = -dtu (8.4) 



where 

3n{t) = J dW")*(x).js(x,t) (8.5) 

is the inner product of the source current with the mode 
eigenf unction. For a compact transmitter loop Ct with 
Nt windings carrying current Ixit) (see Fig. (i for an 
example), this reduces to the line integral [18] 

jn{t) = Nranloit) 

an = (f e(">(x)-dl (8.6) 

Jct 

The solution to ( ) is 

Anit) = Nranlnit) 

In{t) = - / dt'e-^"(*-*')a*-/o(<') 

J —OQ 

= /„(tp)e-^"(*-'-), (8.7) 

where the second line is valid during any quiescent inter- 
val following the termination of a pulse at time tp. Note 
that for a very sharp pulse termination, on a time-scale 
shorter than 1/A„, one may approximate —dt'Io{t') ~ 
AIod{t — tp), where A/q is the down-step size (e.g., about 
5 a in Fig. '.'<). The termination then contributes A/q to 
In{tp)- In any case, it follows that the required coefficient 
in (1.1) is given by 

An = An{tp) = Nranlnitp) (8.8) 

For a compact receiver loop Cr with Nj^ windings, the 
voltage is given by the line integral 

V{t) =NrI E(x,t) -dl. (8.9) 
Jcr 

Substituting (8.3), one obtains the voltage series (1-2) 
with 

Vn = AnNiihn = NtN Ranbnln{tp) (8.10) 

where 

hn = / e(")(x) -dl. (8.11) 
Jcr 

is the corresponding receiver loop line integral. 

This completes the specification of the measured volt- 
age in terms of modal quantities given the transmitter 
and receiver loop characteristics. Note that computation 
of a„ and 6„ requires evaluation of the external electric 
field. For this purpose, the right hand side of (2.41) is 
evaluated using the previously computed internal forms 
(2.42) for the mode eigenf unctions. This evaluation re- 
quires the various quantities worked out in Sec. V for pos- 
itive values of the parameter A(x). The far field asymp- 
totic forms described in (\TI D) may be used at sufficient 
target standoff. The line integrals defining a„, 6„ are then 
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FIG. 7: (Color online) Comparisons of data (solid lines; taken by the NRL-TEMTADS platform [22, 23]) and ah initio theoretical 
predictions (dashed lines) for a 5 x 5 x 10 cm radius aluminum prolate spheroid (left) and a 10 x 10 x 4 cm radius aluminum 
oblate spheroid (right), at various depths-to-center dc, with various axial tilt angles 9. The target is centered under the 
transmitter-receiver pair. Plotted are the measured receiver voltages beginning immediately after the pulse termination times 
seen in Fig. (3, averaged over about ten cycles. The multipliers indicated in each legend entry reflect a ~ 10% variability in the 
transmitter current, as well as residual coil characterization [23] and target positioning errors, and are applied to the data to 
optimize the fit. Straight dashed lines show the predicted l/\/f early time divergence [7]. This first principles agreement, over 
nearly two decades in both time and voltage, is remarkable. 



performed numerically [14] through evaluations at a dis- 
crete set of points along the loops Ct,Cu [19]. The pulse 
wave forms Io{t) are typically given by a sequence of re- 
laxing exponentials and linear ramps, for which /„(i) may 
be evaluated analytically [20]. 

All of these algorithms have been implemented numer- 
ically to produce the comparisons now described. Given 
the precomputation of the mode properties (which re- 
quires 10-20 minutes for a given target on a standard 
workstation), computation of the voltage series (1.2) is 
found to take only about 1 s. This speed is critical to 
efficient solution inverse problems which underlie, for ex- 
ample, the UXO discrimination problem. For the latter, 
properties of an unknown target are estimated by search- 
ing over different candidate targets to find the one that 
produces the best fitting voltage curve predictions [21]. 



2. Comparison with NRL-TEMTADS measurements 

In Fig. we show comparisons between the theory and 
data taken on artificial aluminum spheroids using the 
Naval Research Labs TEMTADS platform [22]. The plat- 
form consists of a 5 X 5 horizontal array of 25 independent, 
well calibrated, high dynamic range concentric transmit- 
ter and receiver coils. The coils are square (35 x 35 cm, 
with Nx = 35 windings for the transmitters; 25 x 25 cm, 
with A^fl — 16 windings, for the receivers) and laid out 



fiat in the same plane with 40 cm between centers in each 
direction [23]. The pulse waveform is described in Fig. fi. 
For illustrative purposes, only the strongest signal, from 
the coil under which the target was centered, is shown in 
Fig. [21]. 

The data quality is seen to be very high, with no vis- 
ible noise or distortion over more than two decades dy- 
namic range of voltage and time. The model predictions 
(dashed lines) are seen to very accurately reproduce the 
data, except at very early time where, as explained ear- 
lier, the absence of the contribution of more complex 
modes in the series (12) with decay times faster than 
about 0.25 ms (A„ > 4000 s~^) cause them to fall below 
the data curves (solid lines). In comparison, the fun- 
damental modes here are a few tens of inverse seconds, 
(as can read off Fig. ' for a^/a — 2,0.4 after applying 
the l/aa^ scaling), and hence have decay times compa- 
rable to the extent of the 25 ms measurement window. 
However, as also noted previously, the early time por- 
tion of the data follows the predicted divergence 
(ultimately cut off only by the finite 10 fis pulse off- 
ramp width, which lies invisibly below the TEMTADS 
measurement window) [7]. Significantly, the N = 7th 
order mean field prediction succeeds in overlapping this 
regime, so that by simply substituting a Vprcd(^*)\/i*A 
tail to the early time prediction (below an optimally cho- 
sen crossover time t* ~ 0.25 ms), one obtains a model fit 
to the that is accurate (at the ~ 5% level) over the entire 
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5 X 5 X 10 cm radius prolate spheroid: Mode voltage arrplitudes vs. tilt angle 
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10 X 10 X 4 cm radius oblate spheroid: Mode voltage amplitudes vs. tilt angle 
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FIG. 8: Plot of the largest individual mode voltage ampli- 
tudes Vn in the series (1.2), plotted vs. mode decay rate A„, 
corresponding to the time-domain curves in Fig. 7 (top: pro- 
late spheroid; bottom: oblate spheroid). Of particular note 
is the complete complementarity of the excited modes for ver- 
tical {6 = 0°; red dots) and horizontal {6 = 90°; blue dots) 
targets, which is primarily responsible for the decay rate of 
the curves at later time. A mixture of the two sets is excited 
at intermediate tilt angles (blue and green circles). It is the 
the failure of the amplitudes to decrease with increasing de- 
cay rate that is responsible for the diverging signals at early 
time. 



dynamic range of the data. 

The different curves in Fig. 7 correspond to a combi- 
nation of different target depths and orientations. The 
changing depth-to-center dc has week effect on the shape 
of the curves, mainly changing the overall amplitude 
(with a roughly l/d^ dipolar dependence). This is a 
reflection of the fact that, at these standoffs and for a 
centered target, the applied magnetic field in the target 
region is fairly uniform and vertical. 

The effect of orientation is more interesting, visibly 
changing the steepness of the curves at later time, with 
the oblate spheroid (right panel) showing a stronger ef- 
fect of this type than the prolate spheroid (left panel). 
The reason for this (as quantified in Fig. , which shows 
a plot of the individual mode amplitudes) is that for a 
vertical symmetry axis {6 = 0), it is horizontally circu- 
lating modes of the type shown in Figs. 4 and ■» that 
are most strongly excited [consistent with the line inte- 
grals (8.6), (8.11)], whereas for a horizontal symmetry 
axis it is the vertically circulating modes, of the type 



shown in Fig. ''). Looking at Fig. 2 (and also the upper 
panel in Fig. S), one sees that for as /a > 1 the lead- 
ing vertical mode (curve containing the lowest red dot) 
decays slightly more slowly than the leading horizontal 
mode (curve containing the lowest blue circle). This ex- 
plains why the red curve {9 = 0°) in the left panel of 
Fig. is slightly steeper at late time than the cyan curve 
{9 — 90°). Conversely, for a^/a < 1, it is the horizontal 
modes that are more slowly decaying (see also the lower 
panel in Fig. s), and the decay rate gap is much larger 
(diverging as a^/a — ^ 0). This explains why, in the right 
panel, the cyan curve is significantly steeper at late time 
than the red curve. 

There is no distinction between the shapes of the 
curves at early time, other than the overall amplitude of 
the divergence. However, the behavior of the am- 

plitudes in Fig. with increasing decay rate explains the 
origin of this divergence in the mode picture. As shown 
in Ref. [7], immediately after pulse termination, the cur- 
rents form a very thin sheet on the target surface. The 
delta-function-like feature requires a superposition of an 
essentially infinite number of modes (cut off only by the 
ultimately finite pulse off-ramp rate), and one indeed sees 
in Fig. 1^ that the mode amplitudes, if anything, are actu- 
ally growing with increasing decay rate. In fact, as can be 
seen explicitly in the exact solution for the sphere, there 
is an infinite subset of modes (corresponding there to a 
fixed subset angular momentum indices I, m that depend 
on the geometry of the background exciting field) which 
have asymptotically constant excitation Vq, and make a 
voltage contribution 



V{t) « Vc 



oo 

P=PO 



Vn 



t^O, (8.12) 



where Aq ~ f / fiaa^ describes the asymptotic behavior of 
the decay rates for large enough p > pq. Equivalently, one 
expects, for a general target with a sufficiently regular 
surface, that there is a density of states Pe(A) ^ [a 
small fraction of the total, which increases as ptot(A) ^ 
\/A], with constant excitation for large A, which indeed 
leads to V{t) ~ / pe(A)e-^*dA - f/v^ for t ^ 0. 

It is observations such as those above, connecting the 
geometry of the decay curves to the geometry of the tar- 
get, that are critical to a workable inversion scheme [21]. 
It is also clear how data from multiple sensors, or from 
different platform positions, which see different effective 
orientations of the target, can greatly aid in such an ef- 
fort. 

As a final comment, the quality of the fits points to 
an interesting implication regarding the accuracy of the 
higher order mode contributions. It is apparent from 
the exact sphere comparison (blue dots) in Fig. that 
the decay rates, and presumably the mode shapes, can 
be trusted quantitatively only, perhaps, for the first few 
dozen modes (recall, also, that this figure shows only the 
first 125 out of 232 computed modes). However, it is clear 
that their summed contribution to the induced voltage 
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measurement is quantitatively extremely accurate. This 
points to the conclusion that the overall effect of a group 
of modes with similar decay rates (in a density of states 
sense) depends only on the part of the mode Hilbert space 
that they cover, not on the detailed partitioning of that 
subspace between individual modes. One can imagine a 
very complex applied field, generated by an intricate set 
of transmitter coils surrounding the entire target that 
is tuned to excite a single high order mode, for which 
the prediction will badly fail. However, for the relatively 
uniform fields of interest here, the conclusion appears 
valid. 



IX. GENERALIZATIONS OF THE THEORY 

In this final section we describe various generalizations 
of the basic theory. In Sec. IX A we consider permeable 
targets, /i ^ /if,- In Sec. IXC we consider simplifica- 
tions in the high contrast limit 3> 1, relevant to 
ferrous targets where /i//ib — O(IO^). The high mag- 
netic contrast limit parallels in many ways that of the 
high conductivity contrast limit (e.g., it enforces a van- 
ishing surface normal n-H of the internal magnetic field). 
However, there are some surprising subtleties in the ex- 
ternal field computation, which is shown to vanish when 
— >■ 0. A computation of the leading 0(/i;,//i) de- 
pendence of n • H is then required, and we show how to 
accomplish this within the Chandrasekhar formalism. In 
Sec. iX t we consider the computation of the freely de- 
caying eigenmodes. Unlike in the nonmagnetic case (Sec. 
II D) where the eigenmodes follow trivially from the di- 
agonalization of the Coulomb integral operator (2.41), in 
the magnetic case the frequency dependence enters the 
operator in a more complicated way, and a sequential 
search must be performed to find the decay rates. In 
Sec. IX D we consider more realistic target geometries, 
including hollow targets and multiple targets. Finally, 
in Sec. IX E we consider the effects of background per- 
meability variations. We have shown that background 
conductivity variations do not impact an induction mea- 
surement, but even a very small background permeability 
has a strong impact. 

All of these generalizations require significantly more 
work to implement, and their applications to experimen- 
tal data will therefore be described elsewhere. 



A. Generalization to inhomogeneous permeability 

For inhomogeneous permeability, (-'.:3) is replaced by, 
1 



/^bV X 



-V X E,, 1 - K^Eb = S 



X 



V X E - K^E = S. 



(9.1) 



Subtracting the first equation from the second, one ob- 
tains 



V X 



1 

Mb 



V X (A - Ab) 



1 



fJ-b 



X H, (9.2) 



in which the magnetic field has been introduced via 
H = B//1 = (iA:/i)"^V x E, and we have again repre- 
sented the electric field in the form (2.1(J), (2.11) in terms 
of a Coulomb gauge vector potentials and the gradient of 
a scalar potential. We now define the background qua- 
sistatic (symmetric) tensor Green function G/i(x, x') by 



V X 



1 

^J■b 



V X Ga(x,x') 
V-Ga(x,x') 



= Pt5(x-x') 

= 0, (9.3) 



in which Pt = S{x - x')l + VV(47r|x - x'|)"i is the 
transverse (divergence-free) projection of the delta func- 
tion. One obtains the formal solution 



A(x)-Ab(x) = / dVG^(x,x') 
1 



) X H(x') 

^J'b , 



[«;2e(x') - «;gEb(x')] 



(9.4) 



The Green function Ga accounts explicitly for variations 
in the background permeability, and in the non-magnetic 
limit ( ) reduces to ( ' i^). Analytic forms for Ga also 
exist, e.g., for horizontally stratified backgrounds. Since 
/i//ib is typically discontinuous at the target boundary it 
is convenient to eliminate the resulting surface term by 
integrating the H term on the right hand side by parts. 
One obtains, 



A(x) - Ab(x) = J Sx' |^G^(x,x') • [aE(x') - abEb(x')] + - l) I^' >< Ga(x,x')] • H(x')| 



(9.5) 



in which the curl operation acts on the second index of restrict the integral to the target volume Vg. For uniform 
Ga- So far, no approximations have been made, but in 
the high contrast limit one may drop the ebEb term and 
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/ib one may replace 



Ga(x,x') ^ 1 



Mb 



47r X — x' 



(9.6) 



The transverse projection terms do not contribute be- 
cause the right hand side of (!•■-), and aU terms derived 
from it in (9.4) and (9.5), are divergence free. With these 
simphfications, (9.5) now reduces to 

,3 ,[m(x')-HH(x') 



V X / (Tx 



47r|x — x' 



1. Coupled integral equations for E and H 

One may now, in principle, use (!) ■'') or (!).7) as the 
basis for a low frequency perturbation theory. However, 
upon substituting for H, the appearance of the curl of E, 
in addition to E itself, on the right hand side is found to 
decrease numerical stability, and it is preferable to use a 
more symmetric approach in which E and H are treated 
on an equal footing. 

One may obtain a second equation coupling the two 
fields by taking the curl of both sides of (;*.5). How- 
ever, the left hand side then produces the combination 
/iH — /ifcHb rather than H — H;,, which turns out to be 
less convenient. In order to obtain the latter, we con- 
struct an alternative to (!>, • ) by formulating the Maxwell 



equations in terms of H instead of E. One may combine 
the Maxwell equations for H and Ht, in the form 

Vx(H-Hb) = -ifc(eE-efcEfc) 
V- [/i,(H-Hfc)] = -V- [(m-M6)H], (9.8) 

in which the source term has been canceled and the right 
hand side of the second equation vanishes outside of the 
target volume Vg- In order to formulate these as an in- 
tegral equation, define the magnetic field tensor Green 
function by 



Gjf(x,x') = — V X Ga(x,x') 

Mb 



(9.9) 



which obeys 



V X Gh(x,x') Pt(5(x-x') 
V- [/ibGH(x,x')] - 0, (9.10) 

and represents a generalization of the Biot-Savart law. 
Define also a background scalar Green function gn sat- 
isfying 



- V • [MbV<?H(x, x')] = <5(x - x'). 



(9.11) 



Together these can be used to construct a formal solution 
to ( ) in the form 



H(x) - Hfc(x) = j d^x' |^Gjf(x,x') • [aE(x') - afcEfc(x')] - {^l - Mb)V[V'.gH(x, x')] • H(x')| , 



(9.12) 



whose structure may be compared to that of ( ). Once 
again, in the high contrast limit one may drop the f7f,Ef, 
term and restrict the integral to Vg. For uniform back- 
ground one obtains 



g//(x,x') = 



1 



47r/i6|x - x'l 



and (?). '-') reduces to 



H(x)-Hb(x) = — Vx / <f 

c Jv, 



,a(x')E(x') 



47r|x — x' 



(9.13) 



(9.14) 



1 

Mb 



VV- / (fix 



,[M( xO-Mb]H(x') 
47r|x-x'| 



whose structure may be compared to (9.7). 

Equations (!).5) and (').l-), or their homogeneous back- 
ground counterparts (9,7) and (9. 14), are the basic results 
of this section. They provide closed integral equations 



for E,H that generalizes the non-magnetic form ( ). 

Defining the coefficients 



Qee = 4:TTikfib/c, Qeh = ik 
Qhh = i/fJ-b, Qhe = 47r/c 



(9.15) 



and the integral equations may be written in the block 
form 



ikA 








H 






+ 



QeE^I QeH^2 



in which the operators 



(tE 
(M-Mb)H 
(9.16) 



3, F(x') 



(Fx 



47r X - 



/C2[F] = V X /Ci[F] 
1C^[¥] = VV-/Ci[F] 



(9.17) 
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represent the basic Coulomb integral operators. These 
may all be shown to be symmetric as well. 



2. Basis function expansion 

The basis function expansion solution to (!l.7) and 
(9.14) requires now separate field expansions [compare 
(2.25)] 



a(x)E(x) = ^aMZf,(x) 

M 

m(x)H(x) = ^&MZf,(x). 



(9.18) 



M 



The electric field basis functions continue to obey 
the divergence free and Neumann boundary conditions 
(2.26). The magnetic field basis functions obey only 
the first condition (except in the limit fi/^J-b oo — see 
below) [24]. For homogeneous ellipsoids, one may con- 
tinue use z|^p [see (2.39)], but simply drop the (1 — x^) 

factor and use zj^^^ = V x [x'+^^'X/m] in place of Z^^^^^ 
[both still rescaled via (2.32)]. 
Define now the coefficients 

at.L = / d3xZf*(x)-Eb(x) 
bt^L - / d'xZf*(x).H,(x), (9.19) 
and the matrix elements 



- 


L 




x).Zf,(x) 

f7(x) 


- 


L 


^3^zr(x).zif,(x) 

/i(x) 


Rlm = 


L 


d^xZf * • 


^ifZf/] 


Slm — 


L 


d^xZf * • 




Ulm = 


L 


cfxZf* ■ 


x;2[(i-Mfc/M)z^] 


Vlm — 


L 


d^xZf; ■ 


^2[Zf(x')]. 



(9.20) 



In terms of these, equation (9.1G) reduces to the super- 
matrix equation [compare (2.27)-(2.30) for the nonmag- 
netic case]: 













it) 



(9.21) 



Qee'R Qeh'^ \ f ^ 
QheV QhhS j I b 



in which O^, O-^, R are self- adjoint. For a homogeneous 
target, S is self adjoint as well, and = (1 — /Xb/^)V. 



The block matrix on the right hand side of (9.21) may 
then be made self adjoint by reexpressing the equations 
in terms of ^/QheO- ~ A*b/A*)a and \/Qeh^- This is im- 
portant for numerical purposes. 

Since the basis functions for ellipsoids remain polyno- 
mials, for the homogeneous target case the Coulomb in- 
tegrals entering (9,20) may all be performed analytically 
using the techniques described in Sees. IV, V and VI. 



B. High magnetic contrast limit 

For ferrous (e.g., steel) targets one typically finds very 
large permeability contrast ^//i6 — O(IO^). Although 
this is far smaller than the O(IO^) conductivity contrast, 
it will often be the case that 1% accuracy is more than 
sufficient, and it is then advantageous to seek simplifica- 
tions in the /i//if, — > oo limit. 

Estimating the terms on the right hand side of (• V), 
one sees that the ratio of the E term to the H term is of 
order 



Ac M ' 



(9.22) 



in which the the curl operation in K,2 is approximated by 
the inverse target size 1/a, and Ac = ATia^a^ is a target 
characteristic decay rate scale. Therefore, for modes with 
decay rates A/ Ac < {^/ the Q eh 1^2 term dominates 
[25]. 

Estimating the terms on the right hand side of (!i I J) 
requires more care. Nominally, the ratio of the E term to 
the H term is also given by ( - '). However, precisely as 
in the high contrast limit for the E field, the nominally 
diverging H term actually forces the boundary normal 
component H • n to scale with the factor ( ), and the 
resulting term [which takes the form of a gradient, pre- 
cisely as does the 1/kI term in (2.9)] serves to cancel the 
boundary normal component arising from the E term. 

In the simultaneous high conductivity and high mag- 
netic contrast limit, (9.16) may therefore written in the 
remarkably symmetric form 



E 
H 



Hfc 



QBif/C2[MH] 

- Qff_EAC2[crE] 



(9.23) 



in which both and are determined by the condi- 
tion that the boundary normal components of their re- 
spective fields vanish. More explicitly, from ( ) one 
identifies 

= -Qmh^ ■lCl[i^i- ^ibm 

^2^,[Mr')/Mb-l]fi(r')-H(r') 



OK 



47r x — x' 



.3^. V'-H(x') 
47r|x - x'l 



(9.24) 



Note that the second term vanishes identically for a ho- 
mogeneous target. The surface integral term clearly di- 
verges with /i/ lib unless the surface normal n(r') •H(r') = 
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0{fib/n) vanishes. This leads to a finite result for V^h 
that appears to depend on the subleading 0{f^b/^J') de- 
pendence of H. However, self-consistently, this term 
must simply act to cancel the leading order surface nor- 
mal component Hf, -I- Qh e^2 [cE] . This leads to a rela- 
tion for that depends only on the leading form for E. 
This derivation is entirely equivalent to the high conduc- 
tivity contrast limit derived in Sec. in which E was 
separated into inductive and gradient parts. 

1. Basis function expansion 



even under inversion, while IC2 is odd. This implies that, 
even for arbitrary ef,/e, /^b//i, if e*^"^(x); h'^"^(x) form an 
eigenmode, so does e'-")(— x), — h("'(— x). The same sym- 
metry implies that one can always choose solutions with 
a definite parity [27], and this implies that e^"\ h*^") 
have opposite parity. However, in the present case, where 
the two are multiples of each other, it follows that both 
e("),±h("^ are eigenmodes, i.e., the rjn must come in op- 
positely signed pairs, and hence that the A„ are each 
doubly degenerate (beyond any other degeneracies aris- 
ing, e.g., from rotation invariance for spheroidal targets). 



Since both E and H now obey the same boundary con- 
dition, one solves (!).2.j) using identical divergence free 
basis functions = = Zm- The identity (2.23) for 
crE, and the analogous one for yitH, implies that the gra- 
dients are orthogonal to the Z^v/ and the basis function 
expansion (" ' ) now yields the off-diagonal form 



O 

o 



a 
b 



(9.25) 



Qeh^ 
QheV 



in which we note that = V is now self adjoint. This 
may be reduced to the single equation for a: 

Oa = at + QEHVO-^hb + QbhQhb VQ-iVa. (9.26) 

The solution to the eigenvalue problem, in which the 
background fields vanish, may be obtained by first solving 
the generalized eigenvalue problem for V: 



Va„ = rjnOOLr. 



(9.27) 



from which one identifies the electric and magnetic eigen- 
vectors 



a^j 



b„ = rjnQHEOLn = 



1 



(9.28) 



QehVu 

with the consistency condition 

QehQhetiI = 1- (9.29) 
Using (9.15) this determines the decay rates A„ = zw„ 



An — 



47777,2 



(9.30) 



For ellipsoids, the volume Vg is symmetric under in- 
version X — —X. It is easy to check that fCi is then 



2. External field 

There is a subtle problem that arises when one at- 
tempts to compute the external field in the high mag- 
netic contrast limit: we will show that the inductive part 
of E — Ef, vanishes identically. Specifically, whenever the 
boundary normal component of H vanishes, /C2[mH] be- 
comes a perfect gradient. The Aii[(TE] term in ( ) 
produces an 0(/if,//i) inductive contribution, as will the 
leading correction to H. Therefore, the total inductive 
part of the external field is of order /ifc//i, and its accurate 
computation requires the leading correction to H. How- 
ever, it is actually only the correction to the boundary 
normal component that one requires, and we will see that 
this may be extracted directly from in (9.23). We will 
show that the leading correction to A^ also follows only 
from this normal component. 

To see that the EMI contribution of the H field term 
comes directly from its surface normal component, note 
that one may write 



1 



47r X 



V X An 



(9.31) 



where Amon is the vector potential associated with a 
monopole field Bmon = One choice is [26] 



tan(6l/2) - 1 - cos(0) 1 - 

= — ■ ,n\ i—n^ 



z X X 
47r(|x| +z)\x\ 



sm{6) 47r|x| 



(9.32) 



Other choices are related to this one by a gauge transfor- 
mation. Using this form, one obtains for any vector field 
F 
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47r X — x' 



X F(x') = [V X A„,on(x - x')] X F(x') 

- [F(x') • V]A^on(x - x') - V[Ao(x - x') • F(x')] 

- -V • [F(x')A,no„(x - x')] - V[A„on(x - x') • F(x')], 



(9.33) 



where, in the last hne, we have used VAmon ^VAmon and V • B = 0. The divergence dot product in the first 
term acts on the B index. Inserting this resuh into in (9.1G), one obtains, using F = (/i — /ib)H: 



V X / d"x 



^^^ "Pi ^^^^^"'^ = - / d\W) /.,]A_,(x - r')n(r') • H(r') 

47r|X-X'| Jgy^ 

- V / d^'x'ifii^') - A^b]A„,o„(x - x') • H(x'). 



(9.34) 



The second term is a perfect gradient, and therefore does 
not contribute to an inductive measurement. The first 
term depends only on the magnetic field surface normal, 
and therefore vanishes, as claimed, to leading order in 



3. Leading surface normal component computation 

Specializing, for simplicity, to a homogeneous target, 
the full fields E, H may in principle be obtained through 
a perturbation expansion in g = fJ-b/ fJ-'- 



H(x) = Ho(x)+5Hi(x) 



(9.35) 



and similarly for E. With this definition, the leading 
order inductive part of the external electric field is given 
by 



Eind(x) = Eb,i„d(x) + QBBCr/Ci[Eo](x) 



(9.36) 



QEH^ib I d^r'A„,on(x - r')fi(r') • Hi(r'). 

'dVs 



On the other hand, from ('* -4) one obtains 
$h(x)= f d?r 



, ,n(r')-Hi(r') 



47r X - 



(9.37) 



which establishes a relation between and the leading 
nonzero surface field normal component and also shows 
that $ obeys Laplace's equation 

V^*// = 0, X ^ dVs. (9.38) 

On the other hand, the vanishing of n • Hq on dVg 
leads, from the second line of ( ), to a leading order 
Neumann-type boundary condition 

n(r) • V$H(r) = n(r) • H;,(r) + QHEn{v) ■ t2[<y^o]{v), 

(9.39) 

for r € dVs- 

Together, (!).37) and (f).38) uniquely define <^h- In- 
verting (!).37) for n-Hi and inserting the result in (9.3G) 
then finally produces the desired inductive component of 
the external field. 



4- Solution via Chandrasekhar approach 

The computation defined by Sec. IX B 3 can be con- 
veniently implemented for homogeneous ellipsoids using 
the Chandrasekhar methods of Sec. III. We outline the 
approach here. 

We assume that Eq has already been determined, ap- 
proximated as a polynomial of some degree N, by solving 
the leading order equation (!).2r)). The result, along with 
a polynomial approximation for the background field, can 
then be used to compute the right hand side of (9.39) 



/^(x) = QHBan(x)-Vx / 

Jv, 



(fx' 



Eo(x') 
47r|x — x' 



n(x) • Hb(x) 



(9.40) 



in the form of a polynomial of degree iV -I- 2. Here 
n = eaXa/a\ is the un-normalized unit vector. Now, 
since <^h obeys the Laplace equation, it must have a 
spherical harmonic expansion 



$h(x) = ^(j)H,i 



(9.41) 



in which we restrict I < N+2. Solving for the coefficients 
4'H,im requires one to compare the normal derivative of 
(9. 4 I ) to (9 -10) on the boundary. Even though x depends 
on 6, <j) on the boundary of an ellipsoid, the spherical 
harmonic expansion still provides a unique specification 
for a function on dVs, and we therefore need to express 
(9.39) in this form. 

To achieve this, we use the fact that x'Yim is a poly- 
nomial of degree / (see Sec. IX A 2), and defining — 
^^(xq/oc)^ {— 1 on the boundary), we can more gen- 
erally express 



ex'Yi 



Im \ 



(9.42) 
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as a polynomial with some known set of coefficients 
3^Zmp;k- Here, the sum is restricted by the condition 
^1+^2 + ^3 = ^ + 2p. Since the spherical harmonics 
are complete, this relationship is invertible: 



n, then 



ki ko 

1 tXjr\ Jij 



1-2-3 - 51 [^"'ikl'-nP^'^^'^'^C^^)' (9-43) 



where the sum is restricted by the same condition. 

Inserting (9.42) (for p = 0) into ('> -H) and taking the 
normal derivative one obtains 



n(x) • V0h(x) 



Substituting (9.43) into (9.44) and setting ^ = 1 the 
boundary condition takes the form 



Im-l'm' H.I'm' i 



(9.45) 



where = Y.p fhnp and similarly for Hb,i„i and 0//,;™, 
and 



k V a ;,m 

(9.46) 

The sum over k is restricted by fci + fc2 + fcs = /', and 
this also explains the value of the last index on y^^. It 
follows that /' — Z > is even, and the matrix $ has 
an upper triangular-like structure. Therefore, if /(x) is 
a polynomial of degree iV + 2, then fim is nonzero only 
for / < + 2, and it follows as well that one requires 
only I' < N + 2: (f>H,i'm' is nonzero only for I < N + 2, 
and only a finite sub-block of $ needs to be inverted. In 
fact, the upper triangular structure allows the inversion 
of this sub-block to be reduced to a sequence of inversions 
of {21 -I- 1) X {21 + 1) blocks in descending order, I = 
N + 2,N + 1,N,... ,3,2,1. 

The next step is to use (9.37) to determine n(r) •Hi(r) 
in terms of . Expanding 



n(r).Hi(r) = ^/.«r'rz, 



(9.47) 



0o(x) = / dV, 

JdV, I 



where 



n(r')||x-r'| 

= dX=ijd'x'^^^e{^,-e) 

= d^,\p^i^i^(j){^/^i) 
= ^[n + 2-x- V](/.(x) 



(9.48) 



(9.49) 



The last integral is of the type computed in Sec. . Note 
that although (j) has degree n + 2, the operator acting 
on (f) in (9.4n) cancels all terms of this degree, and 4>q is 
actually of degree n. Thus, inserting, via {'■■). the poly- 
nomial representation of (9.47) into (9.37) and applying 
the identity (9.48) one obtains a linear relationship of the 
form 



tlH.lm — ^ 'Hlm-l'm'h^ijmi, 



(9.50) 



in which 'Him;i'm' are known coefficients that follow from 
the matrix y and the results of Sec. \ . Here, both /, I' < 
N + 2, so the relationship can be inverted to obtain the 

'hm- 

Finally, the external field can be computed via (').>!()). 
However, the (|x| -l-z)"^ singularity in Amono means that 
this integral is not in Chandrasekhar form. It is therefore 
better to work from the original form 

Eind(x) =E6and+(9£;£;cr/Ci[Eo]+Q£;HAih^2[Hi], (9.51) 

[where, in an abuse of notation, Ejnd here differs from 
that in (9.3U) by various gradients] in which Hi is any 
field with consistent surface normal values. One choice 
is 



Hi(x) = Y.'^'ch'^i^xl^xl^xl^ 

h[ 



,(1) ^ ^^(1) 



(9.52) 



where h 



(1) 



,(1) 



are the coefficients in the 



equivalent monomial expansion form of (0,47). 



restricting again I < N + 2, the problem then reduces to 
computing the expansion coefficients h\l^ in terms of the 
't'lm- However, the Chandrasekhar method makes this 
very straightforward. Let p(x) be a monomial of degree 



C. Freely decaying mode computation 

The freely decaying modes once again correspond to 
solutions to the homogeneous equations (" ) or (5) 21) 
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FIG. 9: (Color online) Top: Exact analytic results for the de- 
cay rate spectrum of hollow steel spheres (/i — 100, a — 5x 10^ 
S/m) as a function of the void ratio Uh ~ rvoid /''sphorc • The 
decay rates increase rapidly as — > 1. Bottom: Exact 
TDEM voltage predictions using the NRL-TEMTADS plat- 
form for a sequence of hollow steel spheres centered 20 cm 
below the transmitter, including the solid sphere (a^ = 0, 
solid black line). Since the surface geometry of all targets is 
identical, the same early time curve (dashed line) fits all data 
sets. However, the multi-exponential regime begins earlier for 
thinner-shelled targets. 



in which the background terms vanish. Solutions will ex- 
ist only for special values of the frequency w = —iX. For 
nonmagnetic targets (see Sec. II D) the frequency enters 
the eigenvalue equation (2.41) as a trivial multiplier (in 
the coefficient Qee), and the modes are straightforward 
solutions to the generalized eigenvalue equation (2.4.'-)). 

For magnetic targets the frequency enters the Q- 
parameters (9.1">) in a way that cannot simply be fac- 
tored out of the 2x2 block operator. Rather, one obtains 



OA = -H(A)A 



(9.53) 



in which A = (a, b)"^ and O, Ti represent the 2x2 op- 
erators in ( - ). Thus, the modes correspond to spe- 
cial values A„ in which H(A„) has a unit (generalized) 
eigenvalue. To compute these one therefore proceeds as 
follows. First, solve the generalized eigenvalue problem 



OA„(A) = 77„(A)H(A)A„(A) 



(9.54) 



for each fixed A. The |?7„(A)| are increasing functions of 
A, and beginning with small A, there will be sequence of 
increasing values A„ for which ry„(A„) = 1 [28]. These 
are the sought after decay rates, and the corresponding 
A„(A„) represent the mode shapes. This search requires 
on the order of ten repeated diagonalizations to find an 
individual A„, hence thousands to find the entire spec- 
trum. The task is therefore quite numerically intensive 
(though at least the subblock matrices R, S, U, V do 
not need to be recomputed at each step). Spectra re- 
sulting from such a computation, using the ferrous value 
jjL = 100, are shown in Fig. 10. 

There is an interesting feature in Fig. 10 that deserves 
comment. Unlike for = /if, (Fig. 2) where the larger 
decay rate values are found to progressively fill in as the 
basis function order N increases, the /i/Aib = 100 case 
shown in the figure apparently displays two classes of 
modes that are widely separated in decay rate (at least 
near the center of the plot). As N increases the gap is 
expected to gradually fill in (as indicated by the exact 
sphere results). 

The pattern occurs because of the large value of /i//if,, 
and is explained by examining the mode shapes them- 
selves. As discussed in Sec. IX B, as /i/Aib -> oo the 
magnetic field boundary normal vanishes, H • n — 0. 
In fact, this convergence is conditional. For given finite 
/i//ift there will be more slowly decaying modes for which 
H • n = 0(/ib//t) is indeed small (compared to the trans- 
verse component |n x H|), and the decay rate converges 
to a finite value in the limit /iti//i — > 0. However, there 
are also magnetically polarized modes for which H • n is 
not small compared to |n x H| (e.g., for which H is fairly 
uniform inside the target), and these decay much more 
rapidly. The mode shapes of the higher and lower groups 
of modes in Fig. Ill indeed exhibit precisely this differ- 
ence. The interesting non-monotonic behavior of some 
of the mode curves with increasing aspect ratio proba- 
bly originates from similar effects (and would go away at 
high enough order N). 

It should be noted that this division is directly anal- 
ogous to the distinct inductive and polarization/non- 
inductive responses discussed in Sees. IB and VII A 
that occur for large conductivity contrast, a/ai, ^ 1, 
though the division here is much less extreme because 
^i/^ib < ct/ct6 [29]. 

Note, finally, that the reduction/truncation of the 
magnetic field basis functions described in Sec. IX B 1 to 
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FIG. 10: Decay rate spectrum vs. aspect ratio a — a^/a for strongly magnetic/ferrous spheroids (/i — 100). All 232 decay rates 
computed using 232 basis functions each for the electric and magnetic fields with order A'^ = / + 2p < 7 are plotted for 201 
values in the range 0.1 < a < 10. Blue dots at a = 1 show exact analytic results for the sphere. The corresponding results for 
nonmagnetic spheroids are shown in Fig. 2. As described in the text, the large value of ^ effectively divides the modes into two 
sets, one with much higher decay rates. The apparent mode gap near the center of the figure would be filled in with increasing 
A^ (as would smaller gaps that contain unmatched sphere modes). 



those with strictly vanishing surface normal ehminates 
these magnetic polarization modes at the outset. The 
corresponding ofF-diagonal reduction (".25) of the matrix 
equation also restores a simple dependence of the decay 
rates on the spectrum of namely the inverse quadratic 
relation (9.30). The approximate validity of this relation 
can also be used to speed up the previously described A„ 
search for large but finite fj./fj.b- 



D. More general target geometries 

1. Hollow targets 

The theory presented here is valid even if Vg is not 
simply connected, or even connected at all. However, the 
Chandrasekhar theory is specific to solid ellipsoids (Sec. 
V). Unfortunately, many of the most interesting targets, 
especially UXO, are hollow (UXO shell thicknesses are in 
the neighborhood of 10% of the radius). At early time 
currents reside only near the outer surface of the target 
and the hollow portion plays no role [7, 8]. However, at 
intermediate time and beyond, when the currents pene- 



trate throughout the target, the dynamics is very differ- 
ent. In particular, the leading decay rates scale inversely 
with the shell thickness (see the upper panel of Fig. !); 
analogous also to the scaling with small aspect ratio seen 
on the left hand side of Fig. '2). The result is a strong 
deviation from the solid sphere result that begins earlier 
and earlier as the shell thickness decreases (lower panel 
of Fig. 9). Accurate modeling in this regime therefore 
requires extension of the formalism to hollow targets. 

There are main two issues. The first is the generaliza- 
tion of the polynomial basis functions described in Sec. 
lie 2. The rescaling (2.32) from the sphere to the el- 
lipsoid works only if the hollow portion is concentric and 
geometrically similar to the surrounding target, i.e., with 
axes defined by a/j = aha, > < 1- The underlying 
i = 1 sphere basis functions [first line of (2..')9)] then 
remain unchanged, while the i — 2 basis functions are 
replaced by 

zLt = V X [{x^ al)il - x')x'+'PXU (9.55) 
The factor (a;^ — af^){l — x^) enforces vanishing of 

(2) 

the normal component at both boundaries. The 
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concentric-similar assumption is almost certainly ade- 
quate for most UXO. 

The second issue is the evaluation of the Coulomb inte- 
gral matrix elements ( ) or ( ). Given a polynomial 
p(x), the basic Coulomb integral may be written as the 
difference of two solid ellipsoid integrals 

P(x) = /Ci[p](x) (9.56) 

7y,(a) 47r|x-x'| Jv,{8„,) 47r|x-x'| 

Both may be evaluated analytically using the methods of 
Sec. V for arbitrary (and arbitrary center placement). 

The problem is that, since x is external to the Vs{sLh), 
the second integral is not a polynomial, but a sum of poly- 
nomials multiplied by elliptic functions [themselves de- 
pending on the nonzero function A(x) defined by (3.25)]. 
The second volume integral entering the matrix elements, 
namely of ^'(x), or its derivatives, multiplied by a second 
polynomial, then requires a numerical evaluation. This 
is not a major barrier to implementing the more general 
theory, but it does require some extra numerical effort 
that will be part of future work. 

2. Multiple targets 

The formalism may be also be extended, with a similar 
level of effort, to multiple targets. The individual target 
basis functions are unchanged from those of an isolated 
target. The difficulty again comes from the second vol- 
ume integral entering the matrix elements: the mutual 
induction manifests in the integrals such as 

Ri2 = f d'xp^i^) [ d^x'-^i^^, (9.57) 
J 1/2 Jyi 47r|x-x'| 

in which the external field in a second target volume V^, 
generated by a polynomial function pi in the first target 
volume V"/, is integrated against a polynomial function 
P2 in the second target. Once again, even for simple tar- 
gets where this external field can be expressed in terms 
of elliptic functions, the second integral requires in gen- 
eral a numerical evaluation. For well separated targets, 
the mutual induction will be weak, and a local gradient 
expansion about the center of will provide accurate 
analytic results. However, for closely spaced, strongly in- 
teracting targets, which is the physically more interesting 
regime, such an approximation will break down. 

E. Effects of inhomogeneous background 
permeability 

Beyond the general results (!)..")) and (9.12), the re- 
sults presented in this section are specialized to the limit 
of homogeneous background permeability. The problem 
is that, unlike background conductivity variations, which 



do not affect the inductive response as long as the conduc- 
tivity contrast remains high, even small soil permeability 
variations (around unity) will impact the induced voltage 
via the /ifc-dependence of the magnetic Green function 
(9.3). 

The biggest impact is the direct refiection from the 
ground surface. In the absence of an air-ground perme- 
ability contrast the soil is essentially transparent to the 
low frequency signals considered here. However, even for 
very small air-ground contrast (say, 1% or less), the direct 
refiection can mask the signals (which remain essentially 
unchanged, at the 1% level) from smaller and/or more 
deeply buried targets. The ground response is much flat- 
ter in frequency than that of the target, and so the lat- 
ter can still be distinguished (though even this assump- 
tion can break down if the permeability is frequency- 
dependent, e.g., if the soil magnetic impurities have non- 
trivial magnetodynamics). If the background permeabil- 
ity has slow enough spatial variation, the problem can 
also be ameliorated through a target-absent background 
subtraction. If the permeability is known and the ground 
is flat, this subtraction could also be computed theoreti- 
cally (essentially from the corresponding plane interface 
Fresnel coefficient). Again, for 1% level accuracy, the tar- 
get response may be computed assuming homogeneous 
ground permeability, and simply added to the ground 
refiection. 

For larger permeability contrast, and/or highly vari- 
able soils (e.g., strongly magnetic volcanic soils [30]) the 
target response becomes much more difficult to discern. 
Given accurate prior knowledge of fib{^), there is no bar- 
rier in principle to computing an adequate approxima- 
tion to the magnetic Green function (9.3) in the target 
neighborhood (e.g., as a polynomial- modified Coulomb 
singularity; analytic forms are also available for horizon- 
tally stratified media) and using this to solve for the in- 
ternal field mode shapes. Intuitively, one expects ad- 
justments in the target current patterns due to mutual 
induction with neighboring magnetic impurity concen- 
trations. More difficult are the background and external 
fields, which will experience fib over a larger volume and 
hence require a more comprehensive model to compute. 

In conclusion, there appears to be some interesting 
future work to be performed modeling the basic phe- 
nomenology of the effects of magnetically variable back- 
grounds, but in the presence of high, unpredictable vari- 
ability [30], the background subtraction problem has no 
simple solution — the number of unknown parameters is 
too large to constrain using limited survey data confined 
to above-ground measurements. 
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based on the Chandrasekhar approach laid out in Sees. 



III-V [3]. 

[16] The non-inductive part may in principle be computed 
as well using (2.20). For homogeneous target and back- 
ground, where the right hand side of (2.2w) is supported 
on the target boundary, this may be formulated as a 
Laplace-Neumann problem with boundary values pro- 
vided by the internal field. However, as discussed in the 
Introduction, this case is not physically realistic, and at 
is seldom suflicently constrained to enable accurate com- 
putation of <l?oxt. 

[17] This is the same estimate as would used to determine the 
rate of decay of current flowing in circular loop. 

[18] There is an interesting measurement theory sublety in- 
herent in the computation (8.2)-(8,ll). Both the trans- 
mitter and receiver loops are conducting wires whose con- 
ductivity should, in principle, be included in the conduc- 
tivity distribution o-(x). Since a/ah will again be very 
large inside the wires, this will have a very large effect on 
the mode shape e'"' (x) in this region, and one may ques- 
tion the use of the unperturbed mode shape in (^,5) and 
(.S.(i). The decay rate A„ will be perturbed as well, but 
only in proportion to their mutual inductance as com- 
pared to the self inductance of the target, which is very 
small. The paradox is explained by converting (8.()) and 
(■ ' ) to the area integrals of the magnetic flux: e.g., 
an ~ (A„^/c) jj^^ h'"'*(x) • ndA, where n(x) is the lo- 
cal surface normal, and h'"' = {c/X„^) V X e*""' is the 
modal magnetic fleld. Under the reasonable assumption 
that the self inductance of the wires is very small (which 
would part of the design specification), the contribution 
of the immediate wire neighborhood to the area inte- 
gral is small, and it is therefore dominated by the much 
larger area over which the perturbed and unperturbed 
h'"' are essentially identical. It follows that the line inte- 
gral (8.6) is insensitive to the presence of the transmitter 
wire, even though the local values of e'"' (x) are highly 
sensitive. For any type of measurement, there is always a 
calibration that must be performed to relate the device 
output to the local values of the unperturbed fields in the 
absence of the measurement apparatus. In this case the 
procedure is relatively straightforward. 

[19] Another technical issue that must be addressed in any 
comparison with data is that modes are computed in 
the natural frame of the target principle axes, whereas 
measurements are performed in the laboratory frame, 
in which both the measurement apparatus and the tar- 
get may have some nontrivial orientation. The electric 
field computations are most conveniently performed in 
the natural frame, so once a lab frame evaluation point 
x is determined from the placement of transmitter and 
receiver loops [which also specifies the lab frame orienta- 
tion of line integral increments dl in (8.6) and (8.11)], 
the lab frame electric field is expressed as E'''''(x) = 
RE"''*(R'^x). Here R is the usual 3x3 rotation ma- 
trix that converts the components of a vector x' in the 
natural frame to those in the lab frame via x — Rx'. 

[20] For larger targets, 1/A„ for lower order modes be much 
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larger (100 ms or more) than the interval between peri- 
odic pulses (typically tens of ms), and the integral ('^.7) 
will encompass several pulses. A typical measurement av- 
erages over many pulses, and the resulting I„{tp) is a 
geometric series that may still be summed analytically. 
Applications of the present theory to target inversion 
problems, where, for example, multi-sensor data provid- 
ing multiple "looks" at the target is extremely useful, 
lies beyond the scope of the present article, and will be 
presented elsewhere. 

The TEMTADS platform, besides its use for scientific 
data collection, is designed to be towed behind a vehi- 
cle (hence the full name, Time-domain Electro-Magnetic 
Towed Array Detection System) and is actively used for 
UXO remediation work in this configuration. See, e.g., 
G. R. Harbaugh, D. A. Steinhurst, D. C. George, J. B. 
Kingdon, D. K. Keiswetter, and T. H. BeU, "EMI Array 
for Cued UXO Discrimination" at http: //serdp-estcp . 
org/. 

The transmitter coil windings are actually spread out 
over 7.8 cm in height. However, modeling them as ID 
ideal loops centered at 3.9 cm above their base suffices 
for targets at reasonable standoff (say, centered 20 cm 
or more below the plane of the loop). The receiver coils 
are compact in height, and lie at the bases of the trans- 
mitter coils, so their centers lie 3.9 cm below the effec- 
tive center of the latter. Due to the rapid ~ decay 
of the signals with target depth, these details turn out 
to be quantitatively important: prior to incorporating 
them, the comparisons in Fig. contained 10-20% sys- 
tematic biases that could not be explained by any known 
measurement uncertainty (e.g., the known ~ 10% pulse 
amplitude variability cited in the caption to Fig. ~) or 
calibration error. 

Since one actually knows that ikfiH = V x E, one could 
use the more restricted set of basis functions = 
V X at the outset. However, it is found that better 
numerical stability is obtained by letting the equations 



themselves impose this identity on the solution. 
[25] One typically finds that for slowest decay rate Ai one has 
Ka ~ TT, which yields Ai ~ Ac/tt^. thus Ac is actually an 
order of magnitude larger than the fundamental decay 
rate. Thus, according to (!).22), for fJ-f/J-b ~ O(IO^), one 
will need A/Ai — O(IO^) for the QeeICi term to domi- 
nate. 

[26] If one integrates — fi • V over some surface S the result 
is the monopole flux through that surface. Using Am- 
pere's law, the same result is obtain by integrating the 
tangential component of A around the contour C enclos- 
ing S. The form ( '. > i ) may then be derived by taking C 
to be the circle at fixed polar angle 6, which then equates 
27r| Ao| sin(^) to the enclosed flux [1 — cos(0)]/2. 

[27] This follows very generally from the formal observation 
that the parity operator commutes with H, irrespective 
of the value of g. Therefore the eigenvectors may always 
be chosen to have definite parity. 

[28] In fact, only half the rjn{X) actually cross unity. The 
other half remain negative for all A and correspond 
to unphysical solutions for which the Maxwell relation 
H = ikfiV X E fails. This explains how the doubling 
of the matrix size, in which ('' ! -) is used together with 
(9 ,5) instead of enforcing the Maxwell relation at the out- 
set by substituting it directly into (9.5), does not lead to 
a doubling of the number of modes. This structure can 
be verified explicitly in the exact solution of the sphere, 
where the rj„ come in equal but oppositely signed pairs. 

[29] It is the existence of these magnetic polarization modes 
for large but finite f^/Hb that also leads to the crossover 
at very early time from the magnetic target t~^^^ charac- 
teristic time-domain voltage power law back to the non- 
magnetic t~^^'^ power law [8]. 

[30] A well known example in the UXO remediation com- 
munity is the former Navy test range on the island of 
Kaho'olawe, Hawaii (see, e.g., http: //en.wikipedia. 
org/wiki/Kahoolawe). 



